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Chapter 1 
Introduction 



Disordered materials have been the subject of numerous experimental and theoretical 
works for a significantly long period of time. There are several forms of disorder 
and they can be classified in the following way. Topological disorder is a structural 
arrangement where the long range order of a crystal is missing. This can be consid- 
ered as the definition of amorphous materials, which are non-crystalline materials. 
Other types of disorder also exist. Substitutional disorder is where the units of the 
elementary cell change randomly with other units in the lattice. Vibrational disorder 
is where the positions of the atoms around their equilibrium positions are random due 
to temperature. Finally, some materials have magnetic disorder. Here the underly- 
ing crystal can posses randomly oriented spins at lattice points. Materials which are 
topologically disordered and posses a randomly oriented spin are called spin glasses. 
They must be not confused with "real" glasses. The latter are prepared by quick 
quenching from the melt. The glass transition appears when the supercooled liquid 
has a sudden change in derivative thermodynamic properties. Henceforth, we will de- 
fine the glasses, which exhibit this glass transition. The glass transition temperature 
is always lower than the melting temperature for crystallines and it partly depends 
on thermal history. 

It is highly important that one is aware of the detailed structure of amorphous ma- 
terials, as this determines the main physical properties. The structure of amorphous 
covalently bonded materials is still however not known in detail and lots of physical 
phenomena are yet to be explained. The amorphous structure cannot be understood 
in terms of a unit cell consisting of a few atoms. The lack of periodicity in the mate- 
rial does not allow for the use of elegant theorems as used for crystalline structures. 
The complex nature of the systems also cannot be described by the means of theo- 
ries used for the periodic structures of solids. In amorphous materials the structural 
information can only be determined in a statistical sense. Diffraction experiments 
are employed to retrieve information about atomic arrangements, which are to yield 
the radial distribution functions containing integrated information about the amor- 
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phous structures. Continuous random network models were proposed to describe the 
structure of covalent amorphous semiconductors on an atomic scale. In the past two 
decades, however, vast advances in information technology has allowed scientists to 
use computers in a more effective way to study disordered materials on the atomic 
scale. 

Transport and energy relaxation of charge carriers in disordered solids has been 
the subject of intensive experimental and theoretical research for the past few decades. 
Detailed knowledge of the real underlying structure is not necessary in order to under- 
stand the electronic properties. Instead, simple models (extended Hamiltonians) were 
proposed, which are able to give reasonable descriptions for these physical properties. 
However, the model parameters used in these Hamiltonians could only be determined 
by some minimal information about the structure. It is possible to describe conductiv- 
ity, mobility, thermopower, and other electronical properties of disordered materials 
with these models. 

This thesis focuses on two physical properties of disordered materials: (i) on the 
structure of amorphous carbon and amorphous silicon, and (ii) on the mobility of the 
quasi one-dimensional disordered organic solids. In chapter |^, the description of the 
molecular dynamics simulation growth of amorphous carbon with low energy carbon 
bombardment is given. The atomistic simulation of the bombardment process during 
the bias enhanced nucleation phase of diamond chemical vapor deposition is presented 
in chapter ^. Chapter ^ deals with the study of the growth of amorphous silicon 
structures by the same preparation method as it was applied for the amorphous carbon 
structures (chapter |^). Finally, in chapter ^, the mobility of quasi one-dimensional 
organic semiconductors will be examined via Monte Carlo simulations. The thesis 
concludes with a summary in English, German and Hungarian. 



Chapter 2 



Molecular dynamics simulation of 
amorphous carbon structures 

Amorphous carbon (a— C) has been the subject of numerous experimental and the- 
oretical works in recent decades. The ability of carbon atoms to form sp^, sp^ and 
sp^ bonding configurations leads to a large variation in structure. Various methods of 
preparation produce a wide range of macroscopic parameters for a— C. Tetrahedrally 
bonded amorphous carbon (ta— C) is found to be optically transparent, extremely 
hard and chemically inert whereas a— C which mostly consists of threefold coordi- 
nated atoms has a small gap in the electronic density of states. In amorphous forms of 
carbon the most elementary parameter is the ratio of fourfold, threefold, and twofold 
coordinated atomic sites which determine electronic and mechanical properties of the 
film. 

In this chapter the growth of amorphous carbon thin film on a [111] diamond 
surface is studied with the tight-binding molecular dynamics technique. Six different 
three-dimensional networks were constructed with periodic boundary conditions in 
two dimensions. Time-dependent non-equilibrium growth was studied and densities, 
radial distribution functions, coordination numbers, bond angle distributions, and 
ring statistics were analyzed. 



2.1 Introduction 

2.1.1 Experimental background 

At the early stage of structure investigations the amorphous carbon was studied using 
diffraction techniques. The obtained results show a wide range of structures. One 
of the first experiments was carried out by a Japanese group on an evaporated a— C 
sample using electron diffraction ||kak60|| . From the measured data they proposed 



a microcrystalline model of amorphous carbon structures, consisting of graphite-like 
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and diamond-like regions. The first neutron diffraction measurement was carried 
out by Mildner and Carpenter 



mi 



il82| 



They estimated a 5% diamond-hke atomic 
configuration in their sample. Li and Lannin published another neutron diffraction 
measurement on a sputtered a— C sample | |i9CI| |. The obtained radial distribution 
functions differed in detail from the theoretical models. A small fraction of sp^ atoms 
were found and the chainlike configurations are also present in the disordered a— C 
network. The absence of a significant fraction of ordered sixfold rings was also shown, 
which was observed in glassy carbon. A diffraction study of a small volume of amor- 
phous carbon materials prepared by filtered vacuum-arc method proposed nearly 
100% of the sp^ atomic configuration ||gas91|| . Another neutron diffraction study was 
employed on an anomalous a— C sample |^hi91|| in order to determine the atomic ar- 
rangement and to decide whether the structure has temperature dependence or not 
[kug93|| . The sample contained mostly sp"^ atomic arrangements. One of the most 
recent neutron diffraction measurement was carried out at the Rutherford Appleton 
Laboratory ||gil95|| . About half a gram of tetrahedral amorphous carbon was prepared 
by deposition on glass substrates. 

In the experiments amorphous carbon films are usually prepared by evapora- 
tion growth techniques. During the evaporation, carbon atoms are directed towards 
the target. If they reach the surface they could (i) scatter backwards, (ii) pene- 
trate under the surface atoms, (Hi) start collision-cascade, or (iv) chemically bond 
on the surface, producing the growth of amorphous carbon. The above mentioned 
processes {(i)—(iv)) depend mainly on the environmental conditions. The most im- 
portant parameters are substrate temperature Tgub and the kinetic energy Eheam of 
bombarding carbon atoms. The diamond-like amorphous carbon is usually prepared 
at low pressure (10~^ mbar), with the substrate temperature lower than 80 °C. In 
the experiments, negative bias is usually applied to accelerate the carbon ions from 
the plasma, causing an extreme increase of initial kinetic energy (usually around 
300—500 eV from 5—10 eV) ||milOO|| . However, there are some neutral atoms that 
mainly contribute to the surface growth because their kinetic energy (1—10 eV) is 
too low to penetrate under the surface. 



2.1.2 Proposed models for tetrahedrally bonded amorphous 
carbon 

The diamond-like amorphous carbon, otherwise frequently referred to as tetrahedral 
amorphous carbon (ta— C), has been the subject of several experiments and theoreti- 
cal works in the last decade. This is so, as the material has lots of very good applicable 
parameters. The ta— C itself has a high percentage of sp^ carbon atoms (nearly to 
even 100%) and the density of states has a wide gap, similar to that of a diamond. 
There are two main models which try to describe how the C"*" ions change the sp"^ 
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rich matrix target to sp^ rich film. These two models were proposed empirically from 
the results of experiments and computer simulations. 

The so-called suhplantation model is based on experimental observation and phys- 
ical intuition. It was proposed in 1989 to describe the film deposition technology for 
semiconductors, metals, and ceramics with low energy E = 1 — 1000 eV hyperther- 
mal species [[lif89| , |lif90|| . Later similar models were proposed by Davis ||dav93|| and 



Robertson [ rob93|| . In these models the ions are implanted into the subsurface of 



the film, producing densification and compressive stress in the process. The ion im- 
pact affects processes on three different time scales: (i) a. collision stage involving 
penetration, energy loss and generation of ballistic displacements (~ 10~^^ s), (ii) a 
thermalization stage, in which the excess energy is dissipated and the local structure 
relaxes (~ 10~^^ s), (in) a. long-term relaxation stage in which temperature activated 
processes may occur (> 10^^'^ s). It is postulated that the ions cause local melting 
near the impact site, a phenomenon known as "thermal spike" which anneals the 
surface region and reduces the stress. 



An alternative model was proposed by Marks et al. ||mar96a|| . This suggests that 



the collision induced stress is stochastically localized on the surface and not in the 
subsurface region, arguing that the growth proceeds directly on the surface and not 
in the bulk region. For this, two-dimensional (2D) MD simulations were done on a 
graphite-like substrate with ion energies from 1 eV to 75 eV. It was emphasized that 
the simulations are not an attempt to reproduce the dynamics of a particular mate- 
rial, it should be viewed as a 2D model of a hypothetical material using a realistic 
interatomic potential. It is really hard to accept any result of this kind of simula- 
tions. However, it was found that the stress calculated in this 2D model reaches the 
maximum at 30 eV deposition energy and then gradually decreases at higher ener- 
gies. This transition has not been previously reproduced by any simulation, but is 
commonly observed experimentally. 



2.1.3 Computer simulations 

There is no experimental method for the determination of microscopic structures in 
three dimensions. Efforts to develop simulation techniques for analyzing atomic scale 
structures are therefore continually being made. The research in this area has recently 
focused on the construction of realistic model structures for amorphous carbon using 
computer simulations by help of two basic methods. Both, the Monte Carlo (MC) and 
molecular dynamics (MD) simulations need local potentials that are usually difficult 
to derive yet are essential to simulations. 

The MC method uses a random number generator to discover the phase space. 
For the initial structure a random packed matrix (or continuous random network) 
is usually chosen. Then an algorithm which satisfies the detailed balance is applied. 
The structure having the lowest energy is believed to have the best amorphous struc- 
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ture. A more efficient algorithm was proposed for continuous disordered systems by 
Barkema et al. in 1996 [ |bar96|| , which is called activation-relaxation technique (ART). 



Another simulation method based on MC simulations has also been developed by Mc- 
Greevy and Pusztai | |mcg88|| for the investigation of non- crystalline structures. The 
Reverse Monte Carlo (RMC) scheme is based on the results of diffraction measure- 
ments and this is the only simulation method which is free from the description of the 



interaction between atoms. It was applied for amorphous carbon [ wal98 | providing 
a large scale three-dimensional model. MC simulations are very successful in find- 
ing the minimal energy structures. However, the dynamics of any system cannot be 
followed. 

Molecular dynamics simulation is a technique to compute the equilibrium and 
transport properties of a classical many-body system. MD simulations are in many 
respects very similar to real experiments. When a MD study is done, a sample 
is prepared consisting of N particles and then the Newton equations of motion are 
solved for this system. Time averages for such a system almost represents that of a 
microcanonical ensemble (E,N,V), if one assumes that time average is equivalent to 
ensemble averages. The microcanonical ensemble is inconvenient both for theoret- 
ical analysis and also for comparisons of computed results with experimental data. 
Canonical and other ensembles can be derived using extended Hamiltonians. The 
reader can find more details about MC and MD simulations, for example, in the book 
by Allen and Tildesley ||all90|| , and in the book by Rapaport [|rap95| . 



In the MD simulations two different interatomic potentials are often applied. The 
ffist set of potentials are the classical empirical potentials. The most empirical po- 
tentials fell into two main groups. Pair potentials like Morse or Lennard- Jones 6—12 
interatomic potentials are inapplicable to covalently bonded systems. The second 
group of interactions contain at least one additional term, i.e. three-body potential 
is included | Jter88| , |bre9U|| . This term could stabilize the covalently bonded struc- 



ture. The above two interactions between carbon atoms are the most commonly 
applied classical empirical potentials for large system with several hundred atoms 
]te?88| , ^^au92i ^d93| , jauOOj jlgOq . 



Although simulations based on empirical potentials achieved remarkable success, 
extension of the simulation method to quantum mechanics is important. A more 
accurate approach is to derive the interatomic potential directly from the electronic 
ground state where the total energy functional is calculated using density functional 
theory (DFT) within the local-density approximation (LDA). Unfortunately, DFT 
is very computer time-consuming, it is therefore unthinkable to perform a sepa- 
rate self-consistent electronic minimization at every MD step. The approach called 
Car-Parrinello method ||car85|| overcomes this difficulty and it has been successfully 



applied for MD simulations of covalently bonded amorphous systems yielding valu- 
able structural and electronic information of small systems. Structural simulations 
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were carried out using ab initio or modified Car-Parrinello DFT method for a— C 
^al89| , |mar96b| , |mccOO| . 



Larger systems are necessary to study for understanding complex dynamics pro- 
cesses sucli as growth. An alternative possibility to Car-Parrinello DFT method is to 
find tight-binding (TB) models that could be applied to larger systems (a few hun- 
dred atoms) and could be less time consuming than DFT |[wan89| , |xu92| , ^an89| , |els98|| . 
In the last decade applications of different TB potentials were widespread in the MD 
study of carbon systems ||dra91| , [wan93a| , |wan93b| , |dra94| , |iun94| , [por95| , [k:oh95| , |don98| , 
|uEl98i |cEu99| . 



The above quantum mechanical potentials were efficiently used to simulate the 
bulk of a— C structures. A common feature of the simulations was that each had an 
initial random network which was then annealed and cooled down near room tem- 
perature. Lots of physical parameters were recovered. Every model shows a very 
good agreement in the short-range order measured in these materials. Furthermore, 
other physical properties were compared like electronic density of states and vibra- 
tional density of states. In 1989, a simulation was made by Galli and her colleagues 
|gal89|| . The disordered carbon network was generated by quenching from liquid using 
the Car-Parrinello method. The 54-atom-supercell reproduced well, the short-range 
order of the atomic structure and the electronic structure of a— C. In the work of 
Wang et al. |[wan93a|| the amorphous structure was generated by quenching from liq- 
uid phase. The supercell contained 216 atoms and the ring statistics analysis showed 
that the clusters formed graphite-like sheets at 2.20, 2.44 and 2.69 g/cw? densities as 
well. It was also found that at higher densities the fraction of fourfold coordinated 
atoms increased. In addition, contrary to previous MD quenching [^al89 



twofold 
wan93b 



coordinated atoms were found in the amorphous matrix. Their latter work 
focused on ta— C. The 74% percentage of fourfold coordinated atoms was generated 
in a sample where the initial density was very high 4 g/cur". Note that the density 
of diamond is close to 3.5 g/cw?. At that density the structure consisted 80% of 
sp^ carbon atoms, and then by decreasing the density, the sp^ content decreased to 
74%. Drabold et al. ||dra94|| criticized the transferability of the potential ||xu92|| used 
in the works of Wang et al. | |wan93a| , [wan93 b|| . However, Wang and Ho argued with 
this statement |[wan94|] . Recently, a more accurate form of their carbon potential was 
published ||tan96|| . The above preparation methods do not contain the information of 
the growth of amorphous carbon. 

Atomic scale modeling of ion-beam-induced growth of amorphous carbon was sim- 
ulated by Kaukonen and Nieminen |[kau92| , [kauOOf and by Jager and Albe ||jagOO|| . The 



work of Kaukonen and Nieminen [|kau92|| - redone later ||kauOO|| - report 1—150 eV 



ion beam kinetic energy window simulations. The empirical Tersoff ||ter88|| potential 
described the interatomic interaction between carbon atoms. It was found in the 
latter work that the low energy (10—40 eV) bombardment produces the highest sp^ 
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content in the host matrix, supporting the model proposed by Marks et al. ||mar96a 
However, the sp^ content was half of the experimental observed value. In the work 
of Jager and Albe [|jagOO|| a similar simulation work was carried out with the Ter- 
soff's ||ter89[ and Brenner's | |bre9CI| | potential. They have found that the fraction of 
tetrahedrally coordinated atoms is too low, even if structures with densities close to 
diamond are obtained. However, changing the cutoff distance in Brenner's potential, 
the sp^ fraction close to 82% was found, which is in good agreement with experi- 
ments. This modification has an undesired side effect, it has increased the fraction 
of fivefold coordinated carbon atoms, which was never observed by any quantum me- 
chanical treatment of interatomic potentials. A general problem of classical potentials 
that are fitted to a limited set of materials properties is their accuracy in describing 
processes which cover areas lying far from the fitted points of the potential energy 
hyperface. In the case of a— C (and t— aC) deposition the validity of using classical 
empirical potentials can be argued. 

It means that quantum mechanical potentials are needed to describe interatomic 
potentials. The first pioneer work was done by Uhlmann et al. [ |uhl9(; ]. In this work 
the diamond-like carbon (DLC, ta— C) films were investigated by bombardment with 
neutral carbon atoms a 3 g/cm^ amorphous film. The results supported the subplan- 
tation model. It was shown that for E > 30 eV the deposition is a subplantation 
process, where the carbon atom species penetrate and occupy subsurface positions, 
producing a defective, low density, sp'^ rich surface layer and a dense, sp^-rich layer 
is formed below this surface. Chu-Chun Fu et al. ||chu99|| published a simulation 



method for the film growth over silicon substrate. The silicon surface is assumed to 
enhance hardness through the formation of a large proportion of tetrahedral bonds. 
They published that a mixed thin layer in the surface preferring carbon-silicon bonds 
to carbon-carbon bonds were formed at impinging energy 7.55 eV of carbon atoms. 
This simulation contained only 9 and 16 bombarding carbon atoms. The above two 
results were published very recently. The project which was done in part of this Ph.D. 
thesis was started long before these results were achieved and will be summarized in 
the next section. 



2.2 Modeling the growth of amorphous carbon 

The aim of this chapter is to describe in detail the low energy (1—5 eV) bombardment 
of the target surface. If the carbon atom reaches the surface with such a low kinetic 
energy the probability of its penetration under the surface and the probability that 
it starts a collision-cascade is very low. It can however sometimes scatter backwards 
and can chemically bond on the surface producing the growth of amorphous carbon. 
These statements allow us to study the growth of amorphous carbon in this environ- 
ment. We address here the analysis of the growth of amorphous carbon with quantum 
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mechanical treatment of the interatomic potential on large scale systems with over 
100 atoms. 



2.2.1 Molecular dynamics calculations 

Molecular dynamics simulations were carried out to study the dynamics of the growth 



xu 



was 



process. The tight-binding (TB) Hamiltonian of Xu, Wang, Chan, and Ho 
used to calculate the interatomic potential between carbon atoms. All the parameters 
and functions of these methods were fitted to the results of local density functional 
calculations. In that scheme, the system is described by a Hamiltonian of the form 



Him 



occupied 



2mi 



+ 



n 



(^„|i^TB({?^})|^n)+^.ep({?^}) 



(2.1) 



where {Fj} denote the positions of the atoms (i = 1, 2, N), and pj is the momen- 
tum of the ith atom. The first term in Eq. (|2.1| ) is the kinetic energy of the ions. The 
second term is the electronic band-structure energy calculated from a parameterized 
tight-binding Hamiltonian ifysdrj}). It consists of the sum of the energy eigen- 
values Ei for the occupied part of the electronic band structure, which is calculated 



by the Slater-Koster empirical tight-binding method |^la54|| with the tight-binding 
parameters determined by fitting the calculated electronic structure of carbon. The 
off-diagonal elements of H^b are described by a set of orthogonal sp^ two-center 
hopping parameters, Vssa, Vspa, Vppa, Vpp-K, scaled with interatomic separation r as 
a function s(r) (Eq. p.2|) ; and the on-site elements are the atomic orbital energies 
of the corresponding atom. The third term is a short-ranged repulsive energy. It 
represents the sum of the ion-ion repulsion and the correction to the double counting 
of the electron-electron interaction in the second term. In this scheme, the quantum 
mechanical nature of the covalent bonding among the carbon atoms is explicitly taken 
into account in the electronic structure which is calculated for each of the atomic con- 
figurations in the MD simulations. Details of the tight-binding model are described 



in reference [ |xu92| ]. The interatomic potential describes accurately the energetic, vi- 
brational, and elastic properties of the diamond (fourfold), graphite (threefold), and 
linear-chain (twofold) structures in comparison with self-consistent first -principles 
density functional calculations. The realistic TB potential has already been used for 
preparation of amorphous carbon ||wan93a|| , tetrahedrally bonded amorphous carbon 



wan93b| ] and successfully apphed for preparation of fullerenes ||zha92| , [las98|| 



The orbital basis used is 2s and 2p. The off-diagonal elements for the tight-binding 
Hamiltonian are given in terms of a function of the interatomic distance s(r) 



s(r) 



exp < n 



+ 



(2.2) 
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This function contains four independent parameters, that have been adjusted with 
results of ab initio calculations or experiments. The diagonal elements, as well as the 
hopping integrals for the equilibrium distance for one particular structure have also 
been parameterized, and are given in the tables of reference | |xu92| |. The repulsive 



energy is a function for the sum of pair interactions, of a similar functional form as 
s(r), with again four independent parameters which are also taken from the previous 
reference 





(2.3) 



To use these interactions in a MD simulation it is necessary to introduce a cutoff 
for the interatomic distances. This is particularly important in cases where periodic 
boundary conditions are imposed, so that replicas of the same atom do not interact 
with each other. 

For integration of Newton's equations of motion the velocity Verlet algorithm was 
used with time step of 0.5 fs: 



ri{t + At) = r,(t) + Vi(t)At + ^ai{t)At^ + 0{At^) (2.4) 

v.(t + At) = v.(t) + ^li^lMim^t. (2.5) 

The attractive forces were obtained using the Hellmann-Feynman theorem. The 
Hellmann-Feynman force on atom i is 

occupied 

F, = - ^ (^„|9i/TB({?^})/5?i|^n), (2.6) 
n 

where the sum is over occupied levels. The repulsive forces are directly the gradient 
of Erep- Finally, we would refer to two articles which describe the above tight-binding 



molecular dynamics method in more details: the work of Wang et al. ||wan89|| with 
orthogonal basis, whose scheme is used in this chapter for simulation of growth process 
of a— C and the work of the Paderborn group ||els98| , |E'aOO|| , which uses non-orthogonal 
basis. 



2.2.2 First simulations 

The review of our computer simulations will be given in this section. This project 
started in 1997 and our aim was to find the best treatment of the description of a— C 
growth using quantum mechanical treatment of applied potential. 

First, we applied K substrate which is the simplest method to simulate the 
energy dissipation process during the collision of bombarding atoms and the substrate 
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atoms. The quantity critical distance was introduced, which was chosen 17% larger 
than the first neighbor bond distance in diamond (1.54 A). We used the following 
model. If a bombarding carbon atom arrived in the neighborhood of a substrate atom, 
it became the part of atoms where the kinetic energy was decreased in every time step 
with a given number, which was smaller, but close to one. If the next bombarding 
atom arrived in the neighborhood of the substrate atom or the closeness of the atoms 
whose kinetic energy decreased in every time step, the similar decrease of energy was 
employed. This model can be considered as the cooling of the bombarding atom in 
the region of the impact, decreasing its kinetic energy in every time step. The zero 
temperature of the substrate provides the quickest energy dissipation. The first results 
showed that the low energy (1.29 eV) bombardment gives a more chainlike structure 
and the higher energy bombardment close to 10 eV provides more amorphous network 
structures ||koh98a| , ^cug99| . 



To improve the description of energy dissipation in our computer experiments the 
task of finding the better temperature control of the substrate was addressed. There 
are several opportunities to control the temperature in MD simulation. The Nose- 
Hoover thermostat uses extended Hamiltonian and describes the canonical ensemble 
in the right statistical physics form. In the case where the particle number N in 
the system tends to infinity, it can be shown that the canonical ensemble can be 
realized approximately in equilibrated systems. However, systems with size N ^ 100 
already show very good results. Unfortunately, it turned out that the non-equilibrium 
growth process cannot be simulated with the above description of constant substrate 
temperature | |koh98b| ]. The failure of the use of Nose-Hoover thermostat does not 



mean that the constant temperature cannot be realized in the simulations. 

For the next step we decided to use substrates, with fixed atoms in the bottom 
layer and giving initially K to the remaining atoms in the substrate. The results 
of the simulations provided similar results as the rigid substrate |[kug98|| with 6.5 and 



13 eV bombardment. The snapshots of the final structures are shown on Fig |2.1| . It 
is supported by the results that lower kinetic energy bombardment supports more 
chainlike structures and that larger energy constructs more amorphous structures. 
It was also shown that the kinetic energy of the atoms in the film decreases near 
exponentially after no further bombarding atoms appeared in the film [ [kohOOf . This 
behavior will be discussed in more detail later in this chapter. 

2.2.3 Numerical algorithm and simulation details 

In this section a more realistic model of amorphous carbon deposition will be given 
[[koh01|| . The simulation technique was as follows. An ideal diamond film consisting 



of 120 atoms was employed to model the substrate. The rectangular simulation cell 
was open along the positive 2;-axis ([111] direction in our case) and periodic boundary 
conditions were used in planar, x and y directions. Atoms in the bottom substrate 
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Figure 2.1: The final structures of tlie atom- by-atom deposition. The left model is chain- 
like amorphous structure with deposition energy 6.5 eV. It contains 286 atoms including 120 
atoms in the substrate. The right model shows the structure made after 13 eV deposition. 
It contains 323 atoms. 

layer were fixed at their ideal lattice positions, whereas the rest of the 96 atoms in 
the uppermost substrate layers were allowed to move with full dynamics. To simulate 
the constant substrate temperature the kinetic energy of the moveable atoms in the 
substrate were rescaled in every time step. The velocity Verlet algorithm was used 
to determine the trajectories of the carbon atoms. Time step was chosen to be equal 
to 0.5 fs. 

Before starting the deposition process, substrate was kept at a given temperature 
for 0.5 ps in order to make structural relaxation. The neutral bombarding atoms were 
randomly placed in x and y directions above the substrate as near as possible, but 
not closer to any other atom than the cutoff distance of the potential. The initial 
velocities of bombarding atoms were directed to the substrate with given kinetic 
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energy. There is no thermal equihbrium so initial velocities were chosen using the 
V = y^ ^gfceam — 0.4 X p) slmple relation, where p is a uniformly distributed random 
number between and 1. Directions were determined by ^ = 120° + 60° x p and 
(p = 360° X p, where 6 and are the polar angles and p is again a random number. 

The frequency of the atomic injection was / = 1/125 fs~^ on average: in every 
MD time step a random number was generated and was compared to a given number 
(in the present case it was 0.004), to see if it is greater or lower. In case of the latter, 
a new particle is "created" and it started its motion bombarding the substrate. It 
means that on average in every 0.004~^ = 250th MD step, the random number was 
smaller than 0.004. The time step was equal to 0.5 fs, i.e. 250x0.5 fs = 125 fs was 
the average time period to put a new atom into the atomic current. This flux is 
in the orders of magnitude larger than the deposition rate usually applied in real 
experiments. We found, however, that even at such a high flux the overheating effect 
of the surface (see, e.g., ||k:au92|| ) can be negligible and the low-energy evaporation 



can be modeled. The lower substrate temperatures applied in our simulations results 
in faster energy dissipation, which compensates for the high deposition rate. The 
initial kinetic energies of the atoms in the beam were dissipated by interacting with 
the surface atoms, whereas their trajectories were followed with full dynamics. We 
note here that some atoms were scattered back and needed to be cut from the system. 
This kind of backscattering can be found in the experiments as well. 

2.2.4 Energy dissipation 

In the experiments the magnitude of the atom flux in the atom-by-atom deposition 
techniques is 10^^ cm~^s~^. It means, that on average one atom bombards the ap- 
proximately 10x10 substrate surface in one second. Computers cannot perform 
simulations in the above time scale. The typical time step chosen is around 1 fs 
and the calculations over several picoseconds last several weeks (months) for systems 
having 100—200 atoms and if the interatomic potential is described by quantum me- 
chanical treatment. Due to the problems mentioned above, all the details of the 
growth cannot be realized by the help of computer simulations. However, the short 
time processes, like energy loss and local thermalization can be simulated, as these 
processes are in picosecond range. It is an approximation that the growth process can 
be simulated with the above time step, when the frequency of bombarding is close to 
1/125 fs~^. Kaukonen and Nieminen ||kau92|| reported that an overheating effect can 



sometimes be observed. To decide weather it plays an important role in our quantum 
mechanical treatment of the problem, we performed computer simulations with only 
one bombarding atom over the 80 atoms substrate surface with 100 K. The time step 
At equal to 0.1 fs was chosen. A carbon atom was placed 2.75 A over the substrate. 
The initial kinetic energy of the bombarding carbon atom was chosen to 10 eV and 
has only ^-directional velocity. We performed simulation placing the carbon atom in 
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Figure 2.2: The 2;-component kinetic energy of the bombarding atom. In the figure on 
the left the time interval 0—150 fs, while in the figure on the right the interval 150—300 fs 
can be seen. 



high symmetrical points over the surface and one general place. The Fig. p.2| shows 
the z component kinetic energy changes for the bombarding atom. It can be seen 
that the carbon atom bounced back after 15 fs and loses its initial kinetic energy 
after 40 fs. Then oscillations of kinetic energy around eV with 0.1 eV deviations 
can be found, as it is clear from the Fig. p.2| . This shows that the 1/125 fs~^ flux is 
acceptable in our computer simulations. 



2.2.5 Structural properties and discussion 



Six different structures were made by our method (see section p.2.3|) : models con- 
structed with 25 ps and 40 ps (denoted by L) injection times, and with average ki- 
netic energy Ebeam = 1 eV and 5 eV at Tsub = 100 K substrate temperature (elTlOO, 
e5T100, elTlOOL and e5T100L) and two other models by 25 ps injection, Ebeam = 1 
and 5 eV energy with substrate temperature of 300 K (elTSOO, e5T300). In all cases 
the structures were relaxed for 5 ps after the deposition, i.e. a bombarding atom 
no longer appeared in the system. The final structures of larger models (elTlOOL 
and e5T100L) made by different bombarding energies consist of almost the same 
number of atoms (177 and 172), with a thickness of 12.65 A and 10.98 A. Models 
elTlOO, e5T100, elT300 and e5T300 have 126—129 atoms. Typical central processing 
unit (CPU) times of simulations are about 50—70 days on DEC-Alpha workstations. 
Snapshots of two films (e5T100L and e5T300) are shown in Fig. and |2.4| after 
growth and relaxation. The substrate at Tsub = 100 K remained similar to crystal 
lattice, while atoms left their position and the topological order became broken dur- 
ing growth at Tgub = 300 K. An atom (open circle) in the middle of the amorphous 
part in Fig. ^.4| originally belonged to the substrate. 
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Figure 2.3: A snapshot of the eSTlOOL model is shown after growth and relaxation. The 
substrate (open circles at the bottom) at Tsub = 100 K remained similar to the crystal lattice 
during growth. Black and grey atoms are fourfold and threefold coordinated, respectively. 
The rest of the open circles correspond to twofold and one onefold coordinated atoms. For 
a color version see: ht t p : / / www . phy. bme . hu / ~kohary / aCfigur es . html 
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Figure 2.4: A snapshot of the e5T300 is shown after growth and relaxation. The substrate 
at Tsttfe = 300 K lost even the topological structure of a perfect crystal during growth. An 
atom (open circle) in the middle of the amorphous part in figure originally belonged to the 
substrate. For a color version see: http: / /www.phy.bme.hu/^kohary/aCfigures.htm| 
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Figure 2.5: The temperature versus time during relaxation of elTlOOL model is shown. 
The fit is described by TfUm = (exp(t[/s]/1069.03 + 6.8) + 150) K, which has an incorrect 
asymptotic behavior. 



Time development of relaxations 

After the growth process the films relaxed with full dynamics for 5 ps. During this 
period the substrate was kept at constant temperature while the deposited networks 
were cooling down. The temperature versus time relation of this non-equilibrium 
process is exponential function in the interval [0:5000] fs as shown in Fig. |2]^ for 
elTlOOL model. 

The best fit was done according to TfUmit) = c + exp(a ■ t + b) K. The constant 
c = 150 K was chosen rather than 100 K, which has a false asymptotic behavior, 
but it gave a better fit. Fitting parameters are a = -0.935 ps~^ and b = 6.8671 and 
c = 150 K. This means that the relaxation time r = —l/a= 1069.03 fs ~ 1.07 ps. 

The time dependence of bond length and bond angle deviations were also in- 
vestigated. In crystals, bond lengths and bond angles only have thermal disorder 
components as a result of the thermal fluctuations. Non-crystalline materials have an 
additional broadening contribution to the bond length and bond angle distribution 
resulting from static disorder. First we calculated the standard deviations aitk) at 
discreet times tk in the following way: 



\ 




dAt,) - SlMi)) (2,7) 



where N is the number of bonds and di{tk) is the ith bond distance. Values between 
0.095 A < (Tfe < 0.11 A were obtained for bond length deviations. (The index b 
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Figure 2.6: As shown, an expected cr^'* and ct@ versus time functions have an exponential 
decay. At the beginning of relaxation when the temperature of the film is about 1000 K the 
thermal part of the fluctuation is important while at the end it plays only a minor role. 



is for bond length.) Similar calculations for bond angle deviations provided aQ{tk) 
within an interval of 11.7°-12.8°. The contribution to standard deviation of thermal 
fluctuation was also studied as a function of time using the following term: 



'y"'[tt) 



\ 



N 



(2.8) 



where time average {di) of bond di is used instead of di{T = K). 

As expected and ctq* versus time functions have an exponential decay as 
shown in Fig. |2]^. At the beginning of relaxation, when the temperature of the film 
is about 1000 K, the thermal part of the fluctuation is important. At the end it plays 
only a minor role. The al^ and Uq values as a function of temperature show linear 
relationships, as displayed in Fig. Straight lines do not intersect origin because 
usually (di) 7^ (ij(T = K). 
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Figure 2.7: The a^^ and Cq in function of temperature show linear relationships. Straight 
lines do not intersect origin because usually (di) ^ di{T = 0). 
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Table 2.1: It contains the number of atoms (No.), the percentage of atoms with dif- 
ferent coordination numbers (Z), the average coordination number (Z), the thickness 
Zd, and the density p of different models. The two rows belong to bulk (top) and 
to total (lower) sample, respectively. For bulk models the densities are always larger 
than for the total sample. 



TV /r„ 1 „l 

iVlodei 


iNO. 


'7 O 


Z=3 


Z=4 






p[g/cm-') 


elTlOO 


113 


0.9 


59.3 


39.8 


3.4 


6.5 


3.0 




129 


7.0 


56.6 


34.9 


3.2 


9.5 


2.3 


eSTlOO 


124 


4.8 


58.1 


37.1 


3.3 


7.7 


2.7 




129 


7.0 


55.8 


35.7 


3.3 


10.7 


2.0 


elT300 


121 


5.8 


71.9 


21.5 


3.1 


9.8 


2.0 




130 


7.7 


70.8 


20.0 


3.1 


12.8 


1.6 


e5T300 


92 


0.0 


68.5 


31.5 


3.3 


6.3 


2.3 




126 


11.1 


64.3 


24.6 


3.1 


9.3 


2.1 


elTlOOL 


160 


5.0 


65.0 


30.0 


3.2 


11.0 


2.4 




177 


9.0 


62.1 


27.1 


3.1 


14.0 


2.1 


e5T100L 


151 


2.0 


55.6 


42.4 


3.4 


9.4 


2.7 




172 


4.7 


56.4 


38.4 


3.3 


12.4 


2.3 



Density 

For density calculations, two different volumes were defined: we will refer to them as 
cells for the total sample and for the bulk. The bottom of the cells were 0.77 A lower 
than the bottom carbon atom in the amorphous network. This is almost at half bond 
distance of bonds between the substrate atoms and the bottom carbon atoms in the 
film. The top of the cells were determined differently. For the total sample the top 
of the cell was the highest z-coordinate that occurs for atoms in the network. For 
the bulk, this z coordinate was decreased by 3 A. The x and y size of the cell was 



determined according to two-dimensional periodic boundary conditions. Table ^TT 
contains the densities of different models. Each row is divided into further two rows. 
The upper rows constantly refer to the bulk and the lower rows refer to the total 
sample. For bulk models the densities are always larger than for the total sample. At 
Tsnfe = 100 K bulk densities are between 2.7 and 3.0 g/cvo? except for model elTlOOL 
where the density is equal to 2.4 g/cm^. The reason for lower density is that on the 
top of the network a void appeared, thus drastically decreasing the local density. At 
Tsnfe = 300 K the structures are less dense (2.0—2.3 g/cm^). 

In Fig. |2.8| density profiles of two networks (elTlOO and elT300) are displayed 
perpendicular to the substrate surface. The arrows represent the layer positions in 
a perfect diamond crystal. The difference in temperature between two substrates 
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Figure 2.8: Density profiles of two films, elTlOO (solid line) and elTSOO (dotted line), are 
displayed perpendicular to the substrate surface. Arrows represent the layer positions in a 
perfect diamond crystal. The difference in temperature between the two substrates causes 
a difference between memory effects. 



causes a difference between memory effects. The structure on the surface at 100 K 
has a more pronouncing layering effect than the other network on substrate at room 
temperature. In the first 3—4 A thick layers over the [111] surface the sp^ content is 
high due to the memory effect. In the rest of the bulk the sp"^ content is dominant. 

Radial distribution function 



The two radial distribution functions of model elTlOOL are displayed in Fig. |2.9| . The 
solid line shows the first and second neighbor contributions to RDF before relaxation, 
while the dashed line belongs to these contributions after relaxation. Peaks in the 
final structure seem to be a bit narrower and in the interval of 1.6—2.2 A there are 
less bond lengths after relaxation. Numerical analysis provides the same conclusion. 



In Fig. |2.1U| partial distribution functions Jij of elTlOOL can be seen, where i and 



j denote the coordination numbers of the connected atoms. Table contains the 
average bond distances which are 1.43 A, 1.52 A, and 1.59 A, for J^ ^, J3 4 and 74^4, 
respectively. 

Coordination number 

For bond counting we used a value of 1.85 A as an upper limit of bond length inside 



the first coordination shell. Table |2.1| shows the percentage of different coordinations. 
There is no fivefold atom in the structures. More than half of the atoms have the sp"^ 
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Figure 2.9: Two radial distribution functions of model elTlOOL are displayed here. The 
solid line shows the first and second neighbors contributions to RDF before relaxation, while 
the dashed line shows these contributions after relaxation. First and second neighbor peaks 
in the final structure are a bit narrower. 





Figure 2.10: The total (C-C) and partial (C3-C3, C3-C4 and C4-C4) distribution 
are shown in elTlOOL model. The '3' and '4' denote the coordination number three and 
four. The average bond distances are 1.43 A, 1.52 A, and 1.59 A for C3-C3, C3-C4 and 
C4— C4, respectively. 
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Table 2.2: Average bond distance statistics 



Model 


C— C distance(A) 




C-C 




C3-C3 


<^C3~C3 


C3-C4 


C"C3-C4 


C4-C4 


(JC4:-CA 


elTlOO 


1.51 


0.09 


1.41 


0.04 


1.53 


0.05 


1.60 


0.07 




1.50 


0.10 


1.41 


0.04 


1.53 


0.05 


1.60 


0.07 


eSTlOO 


1.50 


0.09 


1.44 


0.05 


1.53 


0.05 


1.58 


0.06 




1.50 


0.09 


1.44 


0.05 


1.53 


0.05 


1.58 


0.06 


elT300 


1.48 


0.09 


1.43 


0.05 


1.54 


0.06 


1.58 


0.09 




1.48 


0.09 


1.43 


0.05 


1.55 


0.06 


1.58 


0.09 


e5T300 


1.50 


0.09 


1.43 


0.05 


1.54 


0.06 


1.58 


0.06 




1.48 


0.09 


1.43 


0.05 


1.54 


0.06 


1.58 


0.06 


elTlOOL 


1.49 


0.09 


1.43 


0.05 


1.52 


0.05 


1.59 


0.07 




1.49 


0.09 


1.43 


0.05 


1.52 


0.05 


1.59 


0.07 


eSTlOOL 


1.52 


0.09 


1.44 


0.06 


1.54 


0.06 


1.59 


0.06 




1.51 


0.09 


1.44 


0.05 


1.54 


0.06 


1.59 


0.06 



bonding configuration in the models as shown in Table p.l| . The average coordination 
numbers in the bulk were slightly over the graphite coordination. The dominating part 
of fourfold coordinated atoms is near the substrate and the contribution of twofold 
coordinated atoms is due to the atoms at the top of the amorphous network (see Fig. 



2.3 and 2A). The e5T100L model contains a onefold coordinated carbon atom at the 
end of a chain. When we compare the average coordination numbers there are no 
drastic differences among the simulated models. A MD simulation of Kaukonen and 



Nieminen | |kau92|| , using the classical empirical potential of Tersoff, shows the same 
portion of sp"^ but the sp^:sp^:sp^ ratio is slightly different at Tgub = 300 K. At the 
^beam = 1 gV beam energy 5.8:71.9:21.5 is the ratio in the elT300 model, while it is 
19.0:73.2:6.7 in their case. Our sample has a higher portion of diamond-like atoms 
and less twofold coordinated carbon atoms at the same condition. 



Bond angle distribution 

The angles can be analyzed in detail according to the coordination numbers of the 
neighbors and the center (e.g. C3— C4— C3, C4— C4— C4, etc.) From diffraction data 
this causes more difficulties. In Fig. we display the bond angle distribution 



of model elTlOOL, where the central atoms are threefold (C— C3— C) and fourfold 
coordinated (C— C4— C). In graphite-like crystal and in diamond lattice these angles 
are 120° and 109.47°, respectively. 

The average values of partial bond angle distributions and deviations are gathered 
in Table |2.3| and Table |2.4|. Atoms with sp^ local arrangements have nearly the same 
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Figure 2.11: The bond angle distribution of model elTlOOL, where the central atoms are 
threefold (C-C3-C) and fourfold coordinated (C-C4-C). 



average values in the C3— C4— C3, C3— C4— C4, and C4— C4— C4 cases, which are 
around the tetrahedral angle. Clear differences are observed for sp"^ configurations. 
C4— C3— C4 angles have much less average values than the others. The latter average 
bond angles are close to 120°. It means that threefold coordinated central atoms 
with at least two C3 first neighbors are always near planar structures even if the 
third one is C4. We have not investigated dihedral angle distribution because it is 
not so characteristic at a— C than at amorphous silicon or germanium. 



Ring statistics 

In our analysis the definition of the ring is a closed path, which starts from a given 
atom walking only on the first neighbor bonds. Every atom is visited only once. The 
size of the ring is the number of atoms forming the closed path. The ring statistics for 



our models are displayed in Table p.4| . A fourfold ring was found in the sample made 
with -Efceam = 1 eV at Tsub = 100 K substrate temperature but there was no threefold 
ring in our simulation. In the graphitic and diamond network, sixfold, eightfold, etc. 
rings are present. In our models even- numbered as well as odd-membered rings can 



be found. It seems to be clear from the Table 2.4 that lower beam energy results in 



less 4—7 membered rings than at a higher bombarding energy. 
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Table 2.3: The average values of partial bond angle distributions and deviations. 



Model 


C-C3-C 




CC3C 


a 


C3C3C3 


a 


C3C3C4 


a 


C4C3C4 


a 


elTlOO 


118.1 


9.8 


118.3 


9.6 


118.3 


10.2 


113.5 


7.3 




118.1 


9.5 


118.3 


9.3 


118.3 


10.1 


113.5 


7.3 


eSTlOO 


118.0 


10.5 


117.7 


11.1 


119.9 


9.9 


110.8 


8.0 




118.0 


10.5 


117.7 


11.1 


119.9 


9.9 


110.8 


8.0 


elTSOO 


118.3 


9.5 


119.4 


9.4 


116.4 


9.0 


113.9 


11.1 




118.3 


9.6 


119.3 


9.6 


116.4 


9.0 


113.9 


11.1 


e5T300 


118.1 


10.0 


119.5 


9.3 


118.0 


10.4 


109.9 


7.3 




118.2 


9.4 


119.3 


8.9 


117.8 


10.2 


110.0 


7.1 


elTlOOL 


118.6 


9.4 


119.1 


9.2 


118.1 


9.9 


114.9 


9.9 




118.6 


9.4 


119.0 


9.2 


118.1 


9.9 


114.9 


9.9 


eSTlOOL 


117.9 


10.4 


120.2 


10.3 


115.9 


9.1 


112.7 


14.7 




117.8 


10.2 


119.5 


10.2 


115.7 


9.1 


112.7 


14.7 



Model 


C-C4-C 




CC4C 


a 


C3C4C3 


a 


C3C4C4 


a 


C4C4C4 


cr 


elTlOO 


108.9 


11.1 


107.1 


8.8 


110.2 


11.2 


108.1 


11.1 




108.9 


11.1 


107.1 


8.8 


110.2 


11.2 


108.1 


11.1 


eSTlOO 


109.1 


9.7 


110.6 


9.8 


110.3 


10.7 


108.1 


8.8 




109.1 


9.7 


110.6 


9.8 


110.3 


10.7 


108.1 


8.8 


elTSOO 


109.2 


9.7 


109.2 


11.4 


108.3 


9.2 


111.1 


9.1 




109.2 


9.7 


109.2 


11.4 


108.3 


9.2 


111.1 


9.1 


e5T300 


109.0 


10.0 


108.8 


9.5 


109.3 


10.6 


108.9 


9.4 




109.1 


9.8 


109.0 


9.3 


109.1 


10.4 


108.8 


9.3 


elTlOOL 


109.0 


10.6 


107.8 


10.3 


110.3 


10.9 


107.9 


10.6 




109.0 


10.6 


107.8 


10.3 


110.3 


10.9 


107.9 


10.6 


e5T100L 


109.1 


9.8 


105.9 


11.5 


111.3 


9.9 


108.3 


9.3 




109.1 


9.8 


106.2 


11.7 


111.3 


9.9 


108.3 


9.3 
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Table 2.4: The average values of partial bond angle distributions and deviations. 



TV /r„ 1 „l 

Model 


C-C-C 


Ring statistics 




C— C— c 


a 


4 


5 


6 


7 


elTlOO 


112.8 
113.5 


11.5 
12.2 


1 


21 


40 


37 


e5T100 


113.3 
113.6 


11.4 
12.1 





26 


43 


42 


elTSOO 


115.3 
115.6 


11.3 
11.7 





25 


21 


17 


e5T300 


113.8 
115.2 


10.9 
12.2 





18 


32 


20 


elTlOOL 


114.1 
114.6 


11.2 
12.1 


1 


31 


50 


40 


e5T100L 


112.6 
113.0 


10.9 
11.2 





44 


54 


58 



2.2.6 Summary 

In summary, low energy molecular dynamics simulations of atomic beam growth on 
a diamond [111] surface were carried out using two different average carbon beam 
energies Ebeam = 1 and 5 eV and two different substrate temperatures Tsub = 100 
and 300 K during 30 ps. To obtain larger structures the simulations were prolonged 
by 15 ps at Tsub = 100 K. Six networks were prepared by atom deposition with 
periodic boundary conditions in two dimensions and the most important structural 
parameters were compared. Our aim was to construct large scale amorphous carbon 
models of over 100 atoms by quantum-mechanical treatment. The atomic interactions 
were described by a well-tested tight-binding potential. The structures were grown 
by atom-by-atom deposition onto substrates. The growing process was described as 
in real experiments without any artificial model of energy dissipation [ |k:au92| , |kohOO|| , 
however, the deposition rate was much higher than is usually applied in experiments. 
Other studies on the growth process use similar rates. On the basis of our time 
development investigation we expect that much lower deposition rates do not cause a 
significant difference in our models. The influence of our lower substrate temperature 
is that the energy dissipation is faster and this slightly compensates for the high 
deposition rate. Unfortunately, it is clear that real experimental rates cannot be 
implemented in computer simulation in the near future, even if less time-consuming 
classical empirical potentials are employed. 
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Chapter 3 

Atomistic simulation of the 
bombardment process during the 
BEN phase of diamond CVD 



In this chapter the detailed examination of C^^y projectile bombardments onto hy- 
drogenated amorphous carbon (a— C:H) is studied by means of computer simulations. 
This is the first work addressing these kind of processes with a reasonable quantum 
mechanical treatment. We investigated the conditions of hydrocarbon ion bombard- 
ments in an amorphous carbon layer. This involves the penetration depth study of 
different ion species by different bombarding energies. Special attention was paid to 
the critical energy needed for penetration and to the structural rearrangements in the 
amorphous carbon layer caused by the bombardment. We present here the first simu- 
lation with realistic conditions for bias enhanced nucleation microwave plasma assisted 
chemical vapor deposition of diamond used by Katai et al. [[kat99|, |katOOa|, |katOOb| . 



3.1 Introduction 



3.1.1 The bias enhanced nucleation CVD growth of diamond 

Chemical vapor deposition (CVD) techniques are frequently used to prepare different 
(hydrogenated) amorphous carbon films, as well as diamond layers. Bias enhanced 
nucleation (BEN) is used to achieve a higher nucleation density in the latter case. It 
can be assumed, that during the deposition, the different projectile ions and their ki- 
netic energy with which they arrive at the surface, play an important role. The iden- 
tification of the most abundant ionic species and the measurement of their energy 
distributions under typical microwave plasma enhanced CVD diamond nucleation 
conditions are still an important issue to be addressed if one attempts to validate 
any model of the nucleation process. Katai et al. reported a detailed mass selective 
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energy analysis of the incoming ions during BEN nucleation of chemical vapor depo- 
sition of diamond |[kat99| , [k:at004 |katOOb|| . They measured the different ion energies 
and fluxes as a function of bias voltage at 25 mbar and 800 °C. The total flux of 
hydrogen species (about 0.17x10^^ cm~^s~^) was considerably less than for hydro- 
carbons (around 1.1x10^^ cm~^s~^). At zero bias voltage ions have energies between 
5 and 10 eV and the average ion energy increases for most species with increasing 
bias voltage in a linear way in the low bias regime. The overall ion flux, which is 
integrated over the energy distribution and summed for all ion species, increases with 
ascendant bias voltage almost exponentially. However, the fluxes of several species 
show a maximum versus bias voltage between and 200 V. It was shown that about 
85% of the total hydrocarbon flux is determined by C+, CH+, CH^, CH^, CsH^, 
C2Hj^ species, with average kinetic energies between 50 and 70 eV. The full width 
at half of the maximum energy distributions were between 20 and 30 eV. It was 
manifested that the lower limit for nucleation is around 100 V bias voltage, and the 
optimum is around 200 V. It was observed that among CiHy (y = 1, 2, 3) projectiles, 
the flux versus energy of CHg" distribution is broader and shifted to lower energies 
with a total flux only marginally lower than those of 011012 • The 02112" had the high- 
est total flux of all hydrocarbon species. Finally, the above mentioned experiment 
supported the subplantation model proposed by Lifshitz et al. ||lif89|, |lif90| . 



3.1.2 Proposed model for the formation of diamond nuclei in 
the BEN process 

A speculative model to explain the formation of the heteroepitaxial alignment of 
diamond nuclei with the underlying substrate in the BEN process was proposed on 
the basis of results from in situ plasma analytical measurements, the results from 



surface analytical measurements and from literary data |PtatOOb|| . Here we will give 
just a short summary of this model. 

In the initial phase of BEN a polycrystalline /3-SiO 3—4 nm thick layer is formed. 
It is epitaxially oriented with respect to the underlying silicon substrate. This silicon- 
carbon layer is not always observed in the experiments. However, from the point of 
view of the proposed model, its presence or absence is of no importance. In the next 
step of a typical BEN process, an amorphous carbon layer forms on the surface of 
this silicon-carbon layer. The thickness of the amorphous carbon layer is controlled 
by the balance between simultaneous growth and the etching of the hydrocarbon and 
the hydrogen species. The thickness of this layer was found to be around 5—8 A. 

The penetration depth of ions below the surface during the BEN process de- 



pends on their energy. Uhlmann et al. ||uhl97| , |uhl98|| executed molecular dynamics 



simulations to study the penetration of carbon atoms in diamond like amorphous 
carbon films, which are frequently referred to as tetrahedrally bonded amorphous 



3.2. SURFACE MODEL OF HYDROGENATED AMORPHOUS CARBON 



33 



carbon (ta— C). It was shown for carbon atoms having kinetic energies bigger than 
30 eV that the deposition is a subplantation process, where the carbon atom species 
penetrate and occupy subsurface positions. The penetration depth was found to be 
energy dependent. The resuhs of these computer simulations permit the assump- 
tion that bombarding hydrocarbon ion species having a larger kinetic energy than 
30—40 eV in BEN, penetrate 4—8 A below the amorphous carbon surface. Conse- 
quently, hydrocarbon species penetrate into the vicinity of the epitaxially oriented 
silicon-carbon layer (or the silicon surface), but they do not damage its structure, as 
the amorphous carbon layer is sufficiently thin to give protection. However, a large 
fraction of hydrocarbon ions appears close to the interface of silicon-carbon and the 
a— C layer, increasing the local density in that region of the amorphous carbon layer. 

Following this, the formation of diamond nuclei starts. The proposed model for 
the BEN process did not give a qualitative explanation of this mechanism, but two 
main points were emphasized. The diamond nucleation occurs close to the interface 
of SiC layer and a— C layer, and the silicon-carbon layer is protected from ion damage, 
thus it can provide orientation information. After the nuclei is formed, the kinetic en- 
ergies for almost all of the ions are below the displacement energy of diamond (80 eV) 
and they do not damage the nuclei. 

In this chapter, we concentrate on the incubation period of nucleation. We study 
the penetration characteristics of hydrocarbon ion species by means of computer 
simulations. For this, molecular dynamics simulations were performed by a reasonable 
quantum mechanical treatment. First, a realistic surface model of hydrogenated 
amorphous carbon (a— C:H) film is prepared. Then the bombardments of different 
hydrocarbon ions is studied by different energies. This involves the penetration depth 
study of different ion species by different bombarding energies. Special attention was 
paid to the critical energy needed for penetration and to the possible description of 
the formation of diamond nuclei. Finally, we present here the first simulation with 
realistic conditions for BEN as used in experiments. 



3.2 Surface model of hydrogenated amorphous car- 
bon 



3.2.1 Generation of the surface model 

The surface for our molecular dynamics (MD) simulations was generated from a hy- 
drogenated amorphous carbon bulk supercell with an original density of 2.2 g/cm^. 
This structure was prepared by long time annealing of randomly distributed 128 car- 
bon and 69 hydrogen atoms |[fra93| , |fra94| , [iun94|| . The final arrangement of this model 
has 1.6%, 53.1%, 2.3% and 41.4% fraction of sp^, sp"^, sp^^^ and sp^ bonded carbon 
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atoms, respectively.^] The average first neighbor coordination number is 2.88 in the 
cubic cell with an edge equal to 10.66 A. For bond counting we used the value of 
1.85 A as an upper limit of bond length inside the first coordination shell. We fol- 
lowed previous works | ]don98| , |uhl98|| to prepare a surface for C^Hy ion bombardments. 
First, this structure was doubled in one direction, which we refer to as the 2;- axis. 
Then to model the surface, we broke the periodic continuation along the z direction 
to transform the periodically extended cube into an infinite slab with two free sur- 
faces. One labeled as the "top" surface and the other as the "bottom" surface. We 
then cut the sample along the z-axis, 15 A below the top layer, and the atoms below 
this were removed from the sample. In other words, we had a new bottom surface. 
Henceforth, it shall be referred to as the bottom surface. The cut of periodic continu- 
ation leads to dangling bonds on the surfaces. We chose the top surface as the target 
for bombardments. For the determination of the low energy structure of this film we 
undertook the following. The broken a bonds at the bottom were saturated by hy- 
drogen atoms to retain the electronic structure there. To achieve local rearrangement 
on the top surface (similar to the surface reconstruction in crystals), the sample was 
treated under an annealing process at 2000 K over about 1 ps, with time step At 
equal to 41 a.u. (~ 1.0 fs). This equilibrated structure was then subsequently cooled 
down to 30 K within about 1 ps. Then the sample relaxed with conjugent gradient 
method. Finally, this low energy structure was annealed at 1000 K for about 1 ps 
to prepare the substrate for CxHy projectile bombardment. Our purpose was to do 
simulations with realistic conditions for BEN experiments where the temperature is 
around 800 °C. 



3.2.2 Properties of the model 

The final structure of our surface a— C:H film contained 182 carbon and 107 hydrogen 
atoms, 16 of the latter replace the broken a bonds at the bottom. We chose a 11.1 A 
thick moveable part from the top, which contained 190 atoms. It is sufficient to model 
the bombardment experiments because the deepest penetration length is estimated 
to be around 10 A for the projectiles having kinetic energy lower than 100 eV. The 
rest of the slab at the bottom, with a thickness 4 A, were fixed in all simulations. 
The density of the sample was 1.98 g/cm^. The bulk density was equal to pbuik = 
2.05 g/cm^. We will refer to the bulk, as the structure 3 A below the highest z 
coordinated atom and we will refer to the surface, as the structure above the bulk. 
We used periodic boundary conditions in plane, with the original cell size equal to 
10.66 A. The bonding arrangement contained 6.6%, 48.4%, 5.5% and 38.5% fraction 
of sp^, sp^, sp^'^^ and sp^ bonded carbon atoms, respectively. The mass profile in 

'''The introduction of sp^+^ hybrid denotes a carbon atom with a coordination equal to three, but 
the geometry of the bond having a similar structure to the sp^ hybrid with one "dangling bond" in 
the fourth direction. 
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Figure 3.1: The substrate structure: sp^ fraction, and mass density profile. The 
black atoms show carbon atoms with a coordination equal to four (Z = 4). Grey 
atoms have Z = 3, while light grey atoms have lower coordinations. The white 
circles show hydrogen atoms. The color version of the structure can be found at 
bttp : / / www . physik . uni — mar bur g . de/ ^kohary / CxHyb omb . html . The solid line in the fig- 
ure on the right shows the density of the film and the filled diamonds show the sp"^ fraction. 



Fig. 3TT shows that there is not an important deviation from the average density, the 
minimum is about 1.8 g/cm^ at 11 A and the maximum is about 2.4 g/cm^ at 8 A 
below the surface. 



In Fig. |3.2| the total and partial bond length distributions can be seen. The par- 
tial distributions are centered at 1.42, 1.54, 1.57 A, for C3-C3, C3-C4 and C4-C4 
respectively. {Ci—Cj means the distribution of distances between carbon atoms with 
coordination number i and j). The bond length distributions for the first neigh- 
bor distance show wide fluctuations around the crystalline graphite (1.42 A) and 
diamond (1.54 A) bond lengths. 



3.3 CxHy projectile bombardments. The simula- 
tion details 

In our work we aimed to study the ion bombardment in an sp^-rich a— C:II structure 
by means of computer simulations |P^oh_sub|] . We were interested (a) in the determi- 



nation of average penetration depths, (b) in the examination into the possible decay 
of the projectiles, and (c) in the structural change of the target substrate. We studied 
the effects of the bombardment by different bombarding energies and atomic contents 
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Figure 3.2: Total (C— C) and partial (Ci—Cj) bond length distribution for the substrate 
film. (The i and j denote the coordination of the carbon atoms.) 



of projectiles. The experimental results showed that the most important ions during 



the bias enhanced nucleation (BEN) are C 
C2H3 



CH+, CH+ CH+, C,H+, C,H+ and 



^■3 ) 



kat99 



In fact, C2H2 has the highest total flux among hydrocarbon ions, with 
two or three times larger flux compared to any hydrocarbon and carbon ions. 

We report here, the incubation period of the bias enhanced microwave plasma 
assisted CVD (BEN MW— CVD) phase of diamond nucleation growth. We studied the 
ion bombardment with two relevant projectiles: acetylene (C2H2) and methyl (CH3). 
The interaction between the atoms was described by the density functional tight- 
binding (DFTB) method | ^ls98| , |traOU|| , which was effective to describe the ground 



state geometries and physical properties of different carbon systems ||koh95| , |fra95 



and was also applied in the previous carbon bombardment simulations [|uhl97| , [uhl98 
In the work of Uhlmann et al. the effect of neutral carbon atoms on a 



was studied 



C substrate 

in those simulations the conditions were different to 



However 

ours. They simulated the carbon atom bombardment which resulted in creating a 
diamond like amorphous carbon (DLC or ta— C), which has a high percentage of sp^ 
carbon atoms (even near to 100%). We were however interested in finding a realistic 
simulation for BEN MW— CVD nucleation of diamond and checking the validity of 
the model for the description of heteroepitaxial alignment for diamond nuclei in the 
BEN process, as proposed by Katai et al. |[katOOb|| . 

The simulation details are the following. In the first collision stage, which de- 
scribes the penetration and the energy loss of the species, a projectile with given 
kinetic energy was directed towards the substrate. Random numbers were used for 
determination of the coordinates of the mass center and for the local configurations 
of the species. The lowest atom of the projectile was placed consistently 4 A over the 
substrate, where the interatomic potential almost vanishes. Simulations with three 
different kinetic energies were done (20, 40 and 80 eV) for C2H2 and with 20 and 
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40 eV for CH3 to study the energy dependences of ion bombardment. The choice of 
these two projectiles and the kinetic energies were governed by experimental results 
kat99| , IkatOOai [katOOb|| . The 02132" projectile has the largest flux among hydrocarbon 



species. It can be predicted from previous single carbon projectile bombardments that 
30—40 eV kinetic energy will be the critical energy for penetration of different ion 
species under the surface. We have chosen 20, 40 and 80 eV kinetic energies according 
to the above assumptions. The 20 eV energy bombardment describes the effects of 
low energy projectiles, whose energy is much lower than the energy needed for nucle- 
ation processes. The nucleation starts at 100 V bias voltage. The C2H2 projectiles 
have approximately 40 eV and CH3 ion species have around 35 eV kinetic energy in 
this case. At 200 V bias, at the optimal voltage for nucleation, the CH3 projectiles 
still have kinetic energy around 40 eV, but the C2H2 ions species, which possess the 
highest flux amongst all hydrocarbon ions, have kinetic energy close to 80 eV [[kat99|| . 
This energy is the displacement energy in diamond as well. From CiH^ projectiles 
CH3 was chosen because its flux is just slightly lower than those of CiHq_2. Fur- 
thermore, CiH^ projectiles with low hydrogen content (y < 2) will warrant similar 
effects to single carbon bombardment, which was simulated earlier, but at different 
conditions ||uhl98|| . The collision stage was simulated for 7500 a.u. (~ 182 fs). During 



this period the kinetic energy of the substrate was rescaled according to 1000 K by 
a velocity rescaling method which simulates the constant temperature in the sample. 
For the different kinetic energy ion implantations, different time steps were used: At 
equal to 15, 10 and 5 a.u., for Ekm = 20, 40 and 80 eV, respectively. We studied four 
independent orientations of the bombarding species as it is shown in Fig. ^.3| . The 
first orientation of C2H2 is called POINT, as the top view shows the projectile as a 
point. The second orientation of C2H2 is called LINE. In the case of CH3 projectile, 
two different orientations were called PLANE and EDGE, according their orientation 
towards the surface at the start of the bombardments. The reason for choosing these 
orientations is because they are the extreme cases. In experiments projectile orienta- 



tions are a mixture of these extremes. In Fig. |3.4] the time dependence of the kinetic 
energy of the POINT projectile starting with 80 eV can be seen. The carbon atoms 
have approximately 37 eV kinetic energy, which decreases to 1—5 eV after 400 MD 
steps, when the species arrives to the maximal penetration depth in the film. The 
time step for this case was 5 a.u. (~ 0.125 fs). The hydrogen atoms already have low 
kinetic energy at the start (t = 0). 

In the second thermalization stage, the local structure relaxes and the ion energy 
dissipates into the lattice. The projectile atoms become part of the substrate, their 
velocities were also rescaled according to the substrate temperature (1000 K) for 
approximately 0.5 ps (At 0.5 fs). The trajectories of the next projectile were 
started over this final substrate. In summary, ten different projectiles were started 
over the initial substrate for all different orientations and energies. We note here that 




Figure 3.3: The four different projectile orientations. The names are (a) POINT, (b) 
LINE, (c) PLANE, (d) EDGE. We used the same coordination scheme as in Fig. [Q] . For 
a color version see: |http://www.physik.uni— marburg.de/^kohary/CxHybomb.html| 
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Figure 3.4: In the figure on the left, the kinetic energy versus time for the LINE projectile 
starting with 80 eV kinetic energy is shown. The carbon atoms have initiahy approximately 
37 eV kinetic energy, while hydrogens have much lower energy. One MD step is equal to 
5 a.u. in this case. In the figure on the right, the penetration of the atoms of the projectile 
is shown. The surface of the film is at A. 

some projectile atoms were scattered back and some hydrogen atoms were jostled out 
from the film as well. In Fig. |3.3| we demonstrate projectiles over the final (i.e. after 
ten bombarding events) structures for EDGE40, PLANE40, LINE80 and POINT80. 
The central processing unit (CPU) time for the simulations was 12—15 weeks on HP 
ALPHA workstations. 



3.4 Results 

3.4.1 Dissociation of the projectiles 

One of the main interests of our computer experiments was the time evolution study 
of the species. Monitoring the trajectories of the projectiles showed that their dis- 
sociations are the function of their kinetic energy (see Table |3.1| ). For acetylene the 
carbon-carbon bond was almost never broken at low energies (20, 40 eV). That bond 
was dissociated only twice for LINE orientation at 40 eV. At large bombarding ener- 
gies (80 eV), the carbon-carbon bond was broken nine times for LINE orientation and 
only five times for POINT orientation. Yet, in the experiments the C2H2 projectiles 
have random orientations when they arrive at the substrate surface. This implies the 
fact that at high energies (around 80 eV) the carbon-carbon bond in C2H2 is broken 
between 50 and 90%. The number of broken hydrogen-carbon bonds in the projectiles 
increased with the bombarding energy. One can estimate from Table ^Tl] that more 
than 60% of these kinds of bonds are broken in the experiments. It was also evident 
from the simulations, that after the dissociation, the hydrogen atoms like to form H2 
molecules, particularly in the case of CH3. The final configurations of the projectile 
carbon atoms in the film can be seen in Table ^.2| . They mostly have sp^~^^ and sp"^ 
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Table 3.1: The number of broken bonds during the implantation of the projectiles. 
HI (CI) is the deepest hydrogen (carbon) atom in C2H2 species (POINT orientation). 
HI, H2, H3 are the three hydrogen atoms in the CH3 projectile. 



C2H2 


E(eV) 


Hl-Cl 


C1-C2 


C2-H2 




20 


3 





3 


LINE 


40 


6 


2 


8 




80 


8 


9 


8 




20 


1 





3 


POINT 


40 


6 





8 




80 


9 


5 


10 




CH3 


E(eV) 


C-Hl 


C-H2 


C-H3 




20 


1 


7 


5 


EDGE 


40 


8 


9 


9 




20 


7 


3 


5 


PLANE 


40 


7 


7 


8 



hybrids. For methyl species, usually one hydrogen atom stays in a bond with the 
carbon atom. Also, for acetylene, for lower energy bombardment (20 and 40 eV) the 
hydrogen atoms stay bonded with the projectile carbon atom. We conclude from the 
data that when all the bonds are broken in the projectile during the collision phase, 
i.e. when the carbon atoms become naked, they have a lower coordination number in 
contrast to projectile carbon atoms, when not all the bonds are broken. 



3.4.2 Penetration of the projectiles 

It is important to study the penetration of the projectile atoms to retrieve infor- 
mation about processes during nucleations. In our computer simulations 10 A was 
found to be the deepest penetration length in the sample under the surface. This 
demonstrates that our choice of 11 A as the moveable part of the film was reasonable. 
The average of the maximum penetration depths are orientation, energy and projec- 



tile dependent. Table 3^ shows that higher energy provides a deeper penetration 
depth. There is a significant difference between the implantation depth for POINT 
and LINE projectiles, if the average final penetration depth is considered. This dis- 
tance hardly changes for the orientation LINE increasing the subplantation energy 
from 20 to 40 eV. However, in the same conditions, there are significant changes for 
POINT orientation. We declare here, that orientation has a decisive role in the pen- 
etration behavior of the projectiles, however, in experiments the species arrive with 
random orientations to the surface. 
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Table 3.2: The final orientations of the projectile carbon atoms in the films. CI (C2) 
is the deeper (higher) carbon atom in C2H2 species (POINT orientation). 





LINE 


PLANE 




20 eV 


40 eV 


80 eV 


20 eV 


40 eV 




CI 


C2 


CI 


C2 


CI 


C2 


C 


C 




2 


2 


1 


2 


3 


3 


3 


1 




2 


1 


1 




2 


3 




1 






1 


3 


1 


2 


2 


1 


3 


2 

sp 


3 


2 


3 


4 


3 


2 


1 


3 








1 


1 






1 




■? 

sp 


3 


4 


1 


2 






4 


2 






POINT 


EDGE 




20 eV 


40 eV 


80 cV 
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Table 3.3: The average maximal penetration depths for the projectile carbon atoms. 
The average final penetration depths can be seen in the brackets. For the latter case 
only projectile atoms under the original surface were considered. 



C2H2 


E(eV) 


CI (A) 


C2(A) 


Mass Center (A) 




20 


-2.40 ( -1.31) 


-2.64 ( -1.48) 


-2.38 ( -1.41) 


LINE 


40 


-2.83 ( -2.17) 


-2.73 ( -1.98) 


-2.50 ( -2.05) 




80 


-5.83 ( -5.52) 


-4.71 ( -4.96) 


-4.83 ( -3.98) 




20 


-2.71 ( -1.77) 


-2.81 ( -2.31) 


-2.61 ( -2.02) 


POINT 


40 


-4.76 ( -4.80) 


-4.89 ( -4.61) 


-4.51 ( -4.32) 




80 


-6.24 ( -6.09) 


-5.80 ( -4.59) 


-5.59 ( -4.73) 



CH3 


E(eV) 


C(A) 


Mass Center(A) 


EDGE 


20 


-2.84 ( -2.07) 


-2.74 ( -2.13) 


40 


-3.93 ( -2.93) 


-3.30 ( -2.02) 


PLANE 


20 


-2.74 ( -1.83) 


-2.45 ( -1.59) 


40 


-3.43 ( -3.19) 


-2.84 ( -2.17) 



We now turn to the question which was raised by the proposed model for the 
description of diamond nuclei growth in the BEN process ||katOOb|| . In that model 
it was assumed that the hydrocarbon ion species having kinetic energies bigger than 
30—40 eV, penetrate under the amorphous carbon surface in the vicinity of the silicon- 
carbide layer. In the experiments, the thickness of the amorphous carbon layer was 
found to be around 5—8 A, depending on the balance between simultaneous growth 
and the etching of hydrocarbon and hydrogen species. From Table |3.3| it can be 
asserted that for low energy bombardments (with 20 eV), the penetration of the pro- 
jectile ions take place in surface regions. It seems reasonable to affirm that penetration 
in deeper regions starts at around 30—40 eV. For bigger kinetic energies (80 eV), the 
penetration depths of the projectiles are close to 4—8 A, i.e. in the vicinity of the 
silicon-carbon substrate in the experiments. 



3.4.3 Structural rearrangement 

The bombardment of the projectile atoms cause structural rearrangements in the 
substrate and contribute to the formation of diamond nuclei. In the model describing 
the bias enhanced microwave plasma assisted CVD of diamond nuclei, the formation 
of diamond nuclei was not explained because it was beyond of scope of that model 
[[katOObf . It is assumed in the model that the presence of silicon-carbon or silicon 
substrate contribute to the formation of diamond nuclei, as it has orientation infor- 
mation. It is beyond the work of our computer simulation to report on this process 
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here, due to the fact that we did not have a sihcon substrate. In this section we will 
give a possible explanation of the initial process into the formation of diamond nuclei, 
which involves the occurrence of increasing fraction of sp^ atoms in deeper regions 
after the bombardments. 

Taking into consideration the changes in the total film, our results show that the 
thickness of films have grown by 1-3 A in z direction after ten subsequent bombard- 
ments for all different cases. The sp^ content in the films increased with 1—5%, while 
the sp^ content decreased by the same magnitude. There are no significant changes 
in the total density of the films and on average partial first neighbor distances. 

Here, we concentrate on the changes in the different depths of the film. It was 
reported in a similar study by Uhlmann et al. ||uhl98|] that the carbon bombardments 



form an sp^ rich defective surface layer, whose thickness scales with ion energy. Below 
that layer an sp^ rich layer was found. This was one of the main results in their study, 
which supports the subplantation model in the physical vapor deposition method of 
diamond like amorphous carbon films. Nevertheless, the conditions in the BEN pro- 
cess are different. The temperature is around 1000 K compared to the typically low 
temperature (T < 150 K) usually applied during the preparation of diamond-like 
carbon. Furthermore, the plasma contains several different hydrocarbon species in 
BEN, contrary to diamond-like carbon growth, where single monoenergetic particles 
(mostly single carbon ions) are used. However, our computer simulations of bom- 
bardments with C2H2 and CH3 projectiles show similar conclusions to the work of 
Uhlmann et al. ||uhl98| . In Fig. |3.5| the mass density difference and the discrepancy in 



the sp^, sp^ fraction for the final and initial films can be seen for all cases. Structural 
rearrangements even for low bombarding energy with 20 eV can be seen. The main 
contribution is close to the surface in the region between -4.0 A and -2.0 A, measured 
from the top of the initial substrate. The sp^ content increases with 20—30% and 
the sp^ content decreases by the same magnitude after ten bombardment events. For 
40 eV the changes are significant at the deeper region between -6.0 A and -3.0 A. The 
change in the sp^ and sp^ content is broader compared to the change at 20 eV. At 
80 eV, for C2H2 projectile bombardment, the penetration regime moves deeper and 
can be found between -10.0 A and -5.0 A. In all cases the structural rearrangement 
occurs together with the mass density changes. The effect is the largest at 80 eV. 
We note here that the mass density increase over the top of the initial surface comes 
from the slow growth of the top surface. Therefore, changes in Fig. |3]^ over oA are 
irrelevant. 



It is known from the experiments of Katai et al. |[kat99|| that CH3 projectiles 



have 40 eV kinetic energy at 200 V bias, whereas for, C2H2 it is 80 eV. Nevertheless, 
the kinetic energies of the carbon atoms in the projectiles are the same. It is around 
37 eV for CH3 ion species with 40 eV kinetic energy and approximately the same for 
C2H2 projectile. Yet the latter one has two carbon atoms in one projectile. If one 
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Figure 3.5: The depth profiles for amorphous carbon films after ten bombardment events. 
The solid line show the mass density difference profile between the final and the initial films. 
The circles and stars denote the sp^ and sp^ fraction differences, respectively {rj — rjo and 
5 — do). The vertical straight lines show the original top of the initial film. 



3.5. SUMMARY 



45 



compares the effects of these two projectiles at the optimal bias voltage (200 V) for 
the nucleation process, considerable discrepancies are noticeable. The penetration of 
C2H2 projectiles occurs in deeper regions. In addition, the effects of caused structural 
rearrangements are different. The changes are in deeper regions for C2H2 ion species 
in contrast to CH3 projectiles. We can state here that C2H2 projectiles which has the 
largest flux among the hydrocarbon species, strongly contributes to the local density 
increase and structural rearrangement in the 4— 8 A region under the surface. The 
increase of the content of sp^ hybrids at the optimal 200 V bias voltage for diamond 
nucleation in the region 4—8 A under the surface - which is the vicinity of silicon 
substrate and has orientation information in the experiments - can cause the possible 
formation of diamond nuclei in a "percolated way". This suspicion should be verified 
by computer simulation and it is our main goal in the future. 

3.5 Summary 

The structural rearrangements in a— C:H films were investigated under bombardments 
with Cx-Hj, projectiles. Special attention was paid to the study of different projectiles 
at different bombarding energies. Larger kinetic energies produce deeper subplan- 
tations and lead to a higher probability of the dissociation of the projectiles. For 
acetylene the carbon-carbon bond is hardly ever broken at low energies (20, 40 eV) 
and it is broken in 50—90% of cases for larger energies (80 eV). The computer ex- 
periments showed that around 60% of projectile hydrogen-carbon bonds are broken 
during the subplantation. Usually, one of the hydrogen-carbon bond are not broken 
for CH3 and the other two hydrogens form H2 molecules. At the typical 200 V bias 
voltage, which is optimal for diamond nuclei formation in BEN, the C2H2 projectiles 
have 80 eV and CH3 ion species 40 eV kinetic energy. There are significant differences 
in their effects on the structural rearrangements in the films. For C2H2 projectiles, the 
penetration is deeper and these species cause changes in deeper regions, than for CH3 
projectiles. Whereas, the CH3 projectiles do not contribute to the deep changes and 
cause structural rearrangements only in surface regions. The structural changes in 
the film showed that larger kinetic energy ionic species under the average deposition 
depth generate a thicker sp^ rich film whilst the surface mostly consists of atoms with 
sp^ hybrids. The border between these two regimes descend with larger deposition 
energies. The results support the subplantation model. It will definitely be a fasci- 
nating task in the future to do simulations over a silicon substrate and to study the 
effect of this crystalline layer on the formation of diamond nuclei in the amorphous 
film. 
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Chapter 4 



Growth of amorphous silicon. Low 
energy molecular dynamics 
simulation of atomic bombardment 



In this chapter the growth of amorphous sihcon thin films was studied using the 
tight-binding molecular dynamics technique. The atomic deposition program code, 
which had been developed for studying the growth of amorphous carbon structures 
and which was reported in chapter ||, was applied. 



4.1 Introduction 



Amorphous silicon (a— Si) has been the subject of numerous experimental and theo- 
retical works in recent decades |[mor99|| . Several diffraction measurements were em- 
ployed on a— Si samples in order to retrieve information about the atomic arrangement 
[|moss70| , |bar77| , |for89| , |kug89| , |cic90| , ^g93|| , i.e. the radial distribution function (RDF), 
which contains integrated information about the amorphous structures. The results 
of these experiments suggest that the covalently bonded amorphous silicon are not 
completely disordered. The bonds between atoms and the coordination numbers are 
very similar to the crystalline phase. First and second neighbor peaks are broadened, 
but the positions are the same as in the crystalline structure, while the third peak 
disappears in the measured RDF. The absence of the third peak confirms that there 
is no characteristic dihedral angle. 

In 1985, a continuous random network model of amorphous silicon (WWW model) 
was constructed by Wooten et al. ||woo85|| who carried out Monte Carlo (MC) sim- 
ulations using the classical empirical potential of Keating |[kea66|| . This defect-free 
network includes fivefold and sevenfold rings in addition to the sixfold rings of the 
diamond structure. Since then, several computer generated models have been con- 
structed using classical empirical potentials | [kel88| , [lue89| , |ish97| or applying different 
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quantum mechanical methods |^ti91| , |t6t94] , [hen96| , |herOO|| , but still the WWW model 
is considered to be the best three-dimensional atomic scale representation of the 
amorphous silicon network. Dangling bonds or even floating bonds (if they exist) are 
believed to be very rare in this condensed phase of silicon atoms and are not included 
in the WWW model. Recently, an accurate X— ray diffraction measurement ||laa99 



was published in the wide Q range of 0.03—55 A^^. The most important conclusion 
of this measurement, apart from confirming earlier results, is that the coordination 
number is less than four in implanted a— Si sample, which is controversial with the 
WWW model of amorphous silicon. 

An alternative to MC simulations is molecular dynamics (MD). On the basis of the 
Car-Parrinello density functional theory (DFT) method, MD structural simulation 
was carried out for a— Si having 54 atoms ||sti91|| . Larger systems are needed for 
investigating complex dynamics processes, such as, for example those involved in 
growth. A frequently applied method is to find convenient tight-binding (TB) models, 
which could be employed over a hundred atoms and could be less time consuming 
than DFT. In our simulations the TB Hamiltonian of Kwon et al. [ k:wo94 | was applied 
to describe the interaction between silicon atoms. This group developed an excellent 
TB potential for carbon systems and we successfully applied it to the amorphous 
carbon growth |[k:oh01|| . All the parameters and functions of the interatomic potential 
for silicon were fitted to the results of local density functional calculations. This 
TB potential reproduces energies of different cluster structures, elastic constants, the 
formation energies of vacancies and interstitials in crystalline silicon. According to 
the authors' article, the only disadvantage of this TB model is that the bond lengths 
inside small clusters are a bit longer than those derived from ah initio calculations or 
from the experiments. 



4.2 Simulation details 



Rapid cooling of the liquid phase is frequently applied to construct amorphous and 
glassy structures ||sti91| , |hen96| , |ish97| , |yan9^] . The system is cooled down to room 
temperature at a rate of 10^^ - 10^^ K/s. Another possible approach for constructing 
an amorphous structure |P<:ug_t bp| is motivated by the fact that amorphous silicon 
is usually grown by different vapor deposition techniques in laboratories. We have 
developed a MD computer code (see section |2.2.3| ) to simulate the real preparation 
procedure of an amorphous structure, which is grown by atom-by-atom deposition on 
a substrate. In our recent work ||k:oh01|| , the growth of amorphous carbon films were 
simulated by this method. A brief summary of our simulation technique is described 
here (for details, see reference |[koh01|| or section |2.2.3|) . A rectangular diamond lattice 
cell containing 120 silicon atoms was employed to mimic the substrate. There were 
16 fixed atoms at the bottom. The remaining atoms could move with full dynamics. 
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The simulation cell was open along the [111] direction (positive z-axis) and periodic 
boundary conditions were applied in planar x, y directions. The kinetic energy of 
the atoms in the substrate were rescaled at every time step {At = 0.5 fs) in order 
to keep the substrate at a constant temperature. First, the substrate was kept at 
a given temperature for 0.5 ps to enable structural relaxation. In the deposition 
process the frequency of the atomic injection was on average 1/125 fs~^. This flux is 
orders of magnitude higher than the deposition rate commonly applied in experiments. 
However, the low substrate temperatures (100 and 300 K) during simulations cause 
quick energy dissipation and this may compensate the high deposition rate. 

In our flrst set of simulations the initial kinetic energies of the atoms in the above 
case were dissipated by the interaction with surface atoms. To speed up the energy 
dissipation, we have developed a different simulation method as well. The velocities 
of the atoms in the atomic beam were modifled by the following procedure. If a 
bombarding silicon atom arrived closer than the " critical distance" to any substrate 
atom, it became part of the "group" of substrate atoms. In other words, their total 
kinetic energy was rescaled in every time step to a given (substrate) temperature. If 
the next bombarding atom arrived in the neighborhood of any of the atoms in this 
group, the similar procedure was employed. The chosen critical distance was 15% 
larger than the bond length in crystalline silicon (2.36 A). This is a simple model to 
speed up the inelastic collision between bombarding and target atoms and a plausible 
interpretation of the phenomena. Nine different structures have been constructed by 
the methods mentioned above: 

{i) Six models constructed with 25 ps injection time, and with average kinetic 
energy Eheam = 1 eV and 5 eV, at Tsub = 100 K substrate temperature (elT100R5, 
e5T100R5, elTlOORlO, e5T100R10, elT100R20, and e5T100R20). R5, RIO, and R20 
indicate different relaxation times (5, 10, and 20 ps, respectively) for both structures 
after 25 ps injection. 

{a) Three models (elT300, e5T300, elOT300) were constructed at Tsub = 300 K 
substrate temperature using our model for quick energy dissipation. Average kinetic 
energies of bombarding atoms were Efeeam = 1 eV, 5 eV, and 10 eV, respectively. The 
simulation time was 5 ps for these three models. Due to the preparation method the 
temperature of the fllms remained 300 K, therefore no relaxation was applied. 



[Hi) In order to retrieve information on the difference between rapid cooling |sti91 



|ien96| , |ish97|| and atom-by-atom deposition on a substrate (melt-quenching and vapor- 
quenching), we prepared an additional, tenth model (elTlOOQ; Q is for quenching) in 
the following way. The temperature of the film elT100R20 was increased to 3500 K. 
Considering this, as an initial state, the trajectories of the silicon atoms were followed 
by full dynamics for 7.5 ps. The substrate temperature was kept at 100 K again, which 
leads to the cooling of the film above the substrate. After 7.5 ps, the temperature of 
the fllm decreased into the region at around 200 K. This technique can be considered 
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Figure 4.1: The first and second neighbor peaks in the radial distribution functions for 
elTlOOQ model (sohd hue). The dotted hne shows the experimental silicon RDF, reported 
in [kug93|. It is clear that there is not a determined distance, which can be considered as 
the border for the first neighbor shell, and the first and second neighbor peaks overlap. In 
addition, it can be seen that the model potential applied in these simulations overestimates 
the first neighbor distances. The dashed line represents the J4^4 partial radial distribution 
function, which belongs to distances between Z = 4 coordinated silicon atoms, taking the 
first neighbor shell to 2.6 A. 



as the computer simulation of splat cooling, where small droplets of melt are brought 
into contact with the chill-block. 

The main goal of our work was to simulate amorphous silicon structures. 



4.3 Structural properties 

The final structures of different models consist almost of the same number of atoms 
(between 162 and 177), with a thickness of about 13 A (see Table [4.1| ). In Fig. the 
first and second neighbor peaks are shown in elTlOOQ model. Contrary to carbon 
systems (see section ^, fig. |2^) , there is not a clear, determined distance, which can 
be considered as a border for the first neighbor shell in case for a— Si. We have chosen 
2.6, 2.7, 2.82 and 3.0 A as a limit for the first neighbor shell to calculate the statistics 



of the models. Snapshots of films elT100R20 and elTlOOQ are shown in Fig. ^ 
taking the first neighbor shell to 2.6 A. The substrates with Tsub = 100 K remained 
similar to crystal lattice. 
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Table 4.1: The tables below contain the number of silicon atoms (No.) in the films, 
in addition to the percentage of atoms with different coordination numbers (Z), the 
average coordination number {Z) , the thickness and the density p of the films. Two 
rows belong to bulk (top) and to total (lower) sample, respectively. The numbers in 
the top table (Table A) were calculated by taking the first neighbor shell to 2.6 A, 
and in the bottom table (Table B) by taking 2.82 A. 
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142 


7.0 


19.7 


52.1 


19.7 


0.7 


3.9 


9.7 


2.7 




175 


11.4 


24.0 


45.7 


20.4 


0.7 


3.7 


12.7 


2.4 


elTlOOQ 


122 


1.6 


13.9 


74.6 


8.2 


0.0 


3.9 


9.1 


2.4 




166 


4.2 


22.3 


63.3 


8.2 


0.0 


3.6 


12.1 


2.3 



Table B (2.82 A) 


No. 


Z=2 


Z=3 


Z=4 


Z=5 


Z=6 




zd(A) 


p{g/cm^) 


elT100R5 


140 


0.0 


0.7 


55.0 


31.4 


11.4 


4.6 


10.1 


2.5 




166 


0.0 


3.0 


51.2 


39.3 


13.6 


4.6 


13.1 


2.2 


elTlOORlO 


134 


0.0 


1.5 


53.0 


29.9 


11.9 


4.6 


9.8 


2.4 




166 


0.0 


3.6 


47.0 


39.6 


16.4 


4.7 


12.8 


2.2 


elT100R20 


134 


0.0 


1.5 


53.7 


28.4 


15.7 


4.6 


9.8 


2.4 




166 


0.0 


3.6 


50.0 


35.1 


20.1 


4.6 


12.8 


2.2 


e5T100R5 


151 


0.0 


6.0 


31.8 


43.7 


15.9 


4.8 


10.6 


2.3 




162 


0.0 


5.6 


32.7 


47.0 


16.6 


4.8 


13.6 


1.9 


e5T100R10 


152 


0.0 


2.6 


34.9 


40.1 


15.8 


4.9 


10.7 


2.3 




162 


0.0 


3.1 


34.6 


42.8 


17.1 


4.9 


13.7 


1.9 


e5T100R20 


152 


0.0 


2.6 


35.5 


43.4 


15.8 


4.8 


10.7 


2.3 




162 


0.0 


2.5 


35.8 


46.1 


17.1 


4.8 


13.7 


1.9 


elT300 


151 


0.0 


3.3 


32.5 


40.4 


21.9 


4.9 


10.8 


2.4 




177 


0.0 


7.9 


33.3 


45.0 


21.9 


4.7 


13.8 


2.2 


e5T300 


135 


0.0 


0.7 


27.4 


39.3 


23.7 


5.2 


9.4 


2.8 




171 


0.6 


1.2 


30.4 


47.4 


28.9 


5.1 


12.4 


2.5 


elOTSOO 


142 


0.0 


0.7 


28.2 


45.8 


19.7 


5.0 


9.7 


2.7 




175 


0.6 


2.9 


30.9 


52.8 


21.8 


4.9 


12.7 


2.4 


elTlOOQ 


122 


0.0 


0.8 


56.6 


28.7 


12.3 


4.6 


9.1 


2.4 




166 


0.0 


0.6 


54.8 


41.0 


17.2 


4.6 


12.1 


2.3 
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Figure 4.2: Snapshots of elT100R20 (on the left side) and elTlOOQ (on the right side) 
models are shown after growth and relaxation. The substrates (open circles at the bottom) 
with Tgub = 100 K remained similar to the crystal lattice during growth. Light grey, grey 
and black atoms are threefold, fourfold and fivefold-coordinated, respectively. The rest of 
the open circles correspond to twofold- and one one-fold coordinated atoms. For a color 
version see: |http:/ /www.physik.uni— marburg.de/^kohary/aSifigures.htm] 



4.3.1 Time development of relaxations 

After the growth process, when there was no more injected bombarding atoms over 
the network, the films relaxed with full dynamics for 20 ps. The deposited networks, 
which have a temperature gradient in the z direction released a considerable amount 
of heat. This is attributable to the influence of the substrate kept at a constant tem- 
perature. Anomalous relaxation has been observed during quenching. Temperature 
versus time function of this non-equilibrium process shows a stretched-exponential 
form, Tfiim{t) = c + exp {at^ + 6^ K. 

In the case of the elTlOOQ model the fitting parameters are a = -0.46 (t is 
measured in fs) and b = 7.57 and c = 100 K (substrate temperature). The /3 is equal 
to 0.95 (< 1), producing stretched exponential relaxation. For the elT100R20 model, 
if the fit is done in the [0:10] ps interval, we obtain similar fitting parameters; f5 = 
0.88, a = -0.49, b = 7.54, and for e5T100R20 p = 1.00, a = -0.41, b = 7.79. The 



fits are shown in Fig. [4.3| . We note here that the search for the best fit is not a 
well-defined problem in our case. A slight modification in the parameter list can give 
qualitatively almost the same fit. 
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Figure 4.3: The observed stretched exponential functions for temperature relaxations 
in films clTlOOQ (top), elT100R20 (middle) and e5T100R20 (bottom). The best fit was 
retrieved according to TfUmit) = c + exp (^at^ + 6 j K Fitting parameters are a = -0.46, 

-0.49, -0.41 fs-i; b = 7.57, 7.54, 7.79; P = 0.95, 0.88, 1.00 and c = 100 K, for elTlOOQ, 
elT100R20 and e5T100R20, respectively. The /3 < 1 in all cases, producing stretched 
exponential relaxation. 
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4.3.2 Density 

To ignore the effect of the rough surface on the top side of the grown film we identified 
two different cells: bulk and total sample (see section ^ and [[koh01|] ). The bottom 



of the latter was 0.77 A lower than the lowest z coordinate; and the top of the cell 
was the largest z coordinate, among silicon atoms in the amorphous network. For the 
bulk, this definition was slightly modified, by simply decreasing the z coordinate of 
the top of the cell by 3 A. The x and y size of the cell was determined according to 



two-dimensional periodic boundary conditions. Table contains the densities of dif- 
ferent models. Each row is divided further into two rows. The upper ones always refer 
to the bulk and the second rows refer to the total sample. For realistic density calcu- 
lations one should consider only the "bulk densities", which are always larger than for 
the total sample. At Tsub = 100 K bulk densities are between 2.3—2.5 g/cm^, while 
in the quick energy dissipation models {Tsub = 300 K) the structures possess higher 
densities (2.4—2.8 g/cm^). For crystalline silicon the density is equal to 2.33 g/cm^, 
which is lower than the values we obtained for a— Si, i.e. our tight-binding molecu- 
lar dynamics simulation provided more dense structures. There was no significant 
changing in density during relaxations. 



4.3.3 Ring statistics 

A surprising result was found in ring statistics. We defined a ring as a closed path, 
which starts from a given atom walking only on the first neighbor bonds. Every atom 
is visited only once and bonds can exist only between adjacent atoms. The size of 
the ring is the number of atoms forming a closed path. We considered only ring sizes 
smaller or equal to eight but there are not any eightfold rings in the networks. The 
ring statistics are displayed in Table 0|. Networks prepared by our models have a 



significant number of three membered rings, i.e. triangles are present in the atomic 
arrangement. Furthermore, squares can also be found. Most of the theoretical models 
for a— Si do not contain such fractions of the structure. It seems to be an important 
result although a neutron diffraction measurement carried out on a pure evaporated 
amorphous silicon sample, and evaluated by reverse Monte Carlo method had a similar 
conclusion ||kug93|| . Nevertheless, the systematic analysis of the structural data of 



the Cambridge Structural Database suggests that the equilateral triangles and near 
planar squares can also be natural local configurations inside the atomic arrangements 
of a— Si [|k:ug01||. All of the three independent pieces of evidence support the existence 



of triangles in a— Si structures. 

Significant differences can be found in ring statistics between samples constructed 
by rapid quenching and by atom-by-atom deposition. From Table |4.2| it can be seen 
that for the model elTlOOQ (rapid quenching) the number of rings is much lower than 
for elT100R20 and for elTlOOQ, i.e. models prepared by different techniques have 
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Table 4.2: Ring statistics, for different first neiglibor sliells 2.6, 2.7, 2.82 and 3.0 A, 
respectively. 



Model 


Ring statistics (2.59 A) 


Model 


Ring statistics (2.82 A) 




3 


4 


5 


6 


7 




3 


4 


5 


6 


7 


elT100R5 


12 


10 


50 


86 


95 


elT100R5 


63 


35 


115 


138 


202 


el 1 lUUKlU 


25 


10 


61 


85 


93 


el 1 lOORlO 


73 


37 


122 


169 


232 


Ti nnT? on 
ei i iUUrlzU 


9 


14 


68 


89 


102 


Ti nnT? on 
ei i iUUrlzU 


66 


35 


121 


160 


248 


fi^Ti nnR ^ 
co -L luuno 


11 


9 


51 


45 


65 


n^T^ nnR ^ 

eo -L lUUrLO 


62 


53 


151 


163 


247 


aKT^ nnR i n 
eo ± luufLiu 


14 


21 


74 


84 


96 


o^TinnRin 

eo J- luurLiu 


75 


58 


152 


184 


265 


eo 1 iuurtzu 


18 


20 


78 


103 


116 


eo 1 iuurtzu 


81 


36 


128 


162 


229 


ei 1 ouu 


13 


9 


33 


36 


46 


ei 1 ouu 


63 


45 


84 


87 


97 


eo 1 ouu 


13 


8 


44 


31 


33 


r.p;T'?nn 
eo i ouu 


92 


55 


102 


132 


169 


pi nT'^nn 
eiu ± ouu 


7 


8 


16 


19 


12 


eiu 1 ouu 


54 


43 


71 


60 


82 


qi ti nnn 
ei X iuu><^ 


6 


13 


63 


72 


86 


pi TI nnn 

ei 1 iuuv^ 


75 


29 


109 


129 


181 


Model 


Ring statistics (2.7 A) 


Model 


Ring statistics (3.0 A) 




3 


4 


5 


6 


7 




3 


4 


5 


6 


7 


elT100R5 


39 


26 


83 


111 


157 


elT100R5 


80 


45 


120 


174 


235 


elTlOORlO 


53 


20 


91 


134 


155 


elTlOORlO 


91 


50 


135 


181 


273 


elT100R20 


53 


24 


95 


142 


188 


elT100R20 


82 


56 


137 


177 


276 


e5T100R5 


38 


36 


119 


120 


199 


e5T100R5 


95 


84 


169 


230 


323 


e5T100R10 


38 


35 


118 


133 


173 


e5T100R10 


102 


66 


170 


197 


307 


e5T100R20 


53 


28 


107 


130 


165 


e5T100R20 


95 


58 


157 


188 


293 


elTSOO 


37 


26 


61 


67 


84 


elT300 


81 


55 


87 


95 


118 


e5T300 


56 


33 


78 


91 


131 


e5T300 


123 


72 


136 


158 


277 


elOTSOO 


34 


25 


48 


33 


41 


elOT300 


72 


47 


85 


69 


92 


elTlOOQ 


35 


26 


75 


98 


130 


elTlOOQ 


96 


41 


140 


147 


222 



different medium range order. We can also conclude that for the latter model the 
number of rings has decreased compared to the initial elT100R20 model. With the 
exception of the ring statistics which were different, all the other structural parameters 
of these models were similar. 

4.3.4 Coordination number 

It is known from diffractions that the first neighbor coordination number in amor- 
phous silicon is approximately four. We have found that this number can only be 
provided in our models if we choose 2.6 A for the first neighbor shell. Nevertheless, 
the densities of our models are bigger than in crystalline silicon. This means, our 
amorphous models do not have big voids, indicating that the model potential shows 
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more or less properties of amorphous and liquid silicon systems at the same time. 
In the case of taking the first neighbor shell to 2.82 and 3.0 A, we received coordi- 
nation numbers close to five, which is between the average coordination number in 
amorphous {{Z) ^ 4) and liquid {{Z) ^ 6) silicon. 

Figure [4.1| illustrates the first neighbor distances. It is clear that in the case 
of 2.6 A first neighbor shells, the lowest first neighbor distances are dominated by 
Siz=4— Siz=4 distances (dashed line). The distances between other coordinated sil- 
icon atoms only contribute to larger distances than the interatomic separation in 
crystalline silicon (2.36 A) and do not contribute to lower distances, which can be 



seen in the experimental curve (dotted line) in Fig. ^A . 

To our present knowledge, there is no silicon tight-binding model potential, which 
can more accurately describe the amorphous phase. Vink et al. proposed a modified 
Stillinger- Weber classical empirical potential very recently [|vin01|] , which they believe 
can be used to describe the structure of amorphous silicon. 

4.3.5 Radial distribution and bond angle distribution func- 
tions 



In Fig. |4.1| the first and second neighbor peaks of elTlOOQ model can be seen. All 



the other models show similar radial distribution functions. The tight-binding poten- 
tial applied in our computer experiments overestimates the first neighbor distance. 
However, if the potential would provide only silicon atoms with a coordination equal 
to four, the position of the peak for the first neighbor shell would be in the right 
position. 

The main contribution to the bond angle distribution arises from angles between 
95° and 125°. We found a large amount of triangles in our samples, there is therefore 
a significant contribution from bond angles close to 60°. We believe this causes the 
average bond angle to be smaller than 109°, which is the tetrahedral angle in diamond 
structures, in our systems. 



4.4 Summary 

We have developed a tight-binding molecular dynamics computer code to simulate the 
preparation procedure of amorphous silicon, which are grown by a vapor deposition 
technique. Nine structures were simulated by this method. An additional tenth 
model (elTlOOQ) was prepared by rapid cooling in order to make a comparison 
between the atom-by-atom deposition on a substrate and melt-quenching preparation 
techniques. 

The most important difference we have found between the models prepared in 
varied conditions, is in the ring statistics, i.e. in medium range order. In our simula- 
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tions triangles are present in the atomic arrangement. Two other independent pieces 
of evidence support the existence of triangles in a— Si structures. Nevertheless, the 
fragments like equilateral triangles have never been considered at electronic density 
of states calculations or with unsolved problems such as the breaking of weak bonds 
by prolonged illumination, which is the major mechanism for the light-induced defect 
creation |bhi95||. 



The tight-binding model potential used in our computer experiments slightly 
overestimates the bond distances, giving values which were observed previously in 
liquid forms | |ish96| | . A more important problem is that the TB potential provides 
overcoordinated atomic sites (five- even sixfold coordinated atoms) during the amor- 
phous silicon simulations. The sixfold coordinated atoms we obtained can be found 
in the liquid phase of silicon instead of in the amorphous structure. The well-known 
Stillinger- Weber classical empirical potential for amorphous silicon had a similar dif- 
ficulty, yet this has recently been successfully modified |[vin01|| . A possible future aim 
is to improve this TB potential by a similar process. 
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Chapter 5 

One-dimensional hopping in 
disordered organic solids 



5.1 Introduction 



Transport of charge carriers in disordered organic solids, such as molecularly doped 
polymers, conjugated polymers and organic glasses has been the subject of intensive 
experimental and theoretical study for about 20 years ||bas93|| . In recent years partic- 
ular attention has been devoted to the columnar discotic liquid-crystalline glasses due 
to their potential technical applications for electrophotography and as transport ma- 
terials in light-emitting diodes ||bac97| , |chr97a|| . Dielectric measurements have clarified 
that charge transport in most of such materials is extremely anisotropic ||bod95| 



so 



that it is reasonable to describe the transport of charge carriers as a one-dimensional 
process [[ble99|| . Moreover, experimentally observed dependences of the conductivity 
on the frequency, on the strength of the applied electric field and on the temperature 
evidence that an incoherent hopping process is the dominant transport mechanism in 
such materials ||bod95| , [ble99| , |Dch994 |och99b|| . It has been found reasonable to de- 
scribe the transport of charge carriers as a hopping process via molecules arranged in 
spatially ordered one-dimensional chains ||ble99|| . The energy levels which are respon- 
sible for charge transport are usually characterized by a Gaussian density of states 
(DOS) iEas93| , |duH9^ , [BIi99| , |och99a|| : 



V2 



exp 



2a2 



(5.1) 



The origins of the energy disorder characterized by parameter a could be fluctuations 
in the polarization energy as well as dipolar interactions or impurity molecules ||ble99| . 
Such systems represent an excellent field to test various theoretical approaches. These 
approaches differ from each other, for example, with respect to the choice of transition 
rates for single hopping events. The choice of transition rates is of course limited as 
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they must fulfill the detailed balance condition to be physically reasonable in thermal 
equilibrium. The most plausible transition probabilities for the hopping of charge 
carriers via localized states were suggested by Miller and Abrahams [|mil60|, [bas93| , 
f3le99| |: the jump rate from an occupied localized state i to an empty state j separated 
by distance Rij is the product of a prefactor Fi, overlap of the wave functions of states 
i and j, and the Boltzmann factor for jumps upward in energy |[mil60||, 



Fi exp 



2R. 



1,3 



a 



exp 



ej - ei 



2kT 



(5.2) 



Here a is the decay length of the carrier wave function in the localized states and k 
is the Boltzmann constant. Allowing only hops to nearest neighbors and assuming 
that chains of localized states are spatially regular we omit the /^-dependence of the 
transition probabilities and use the Miller- Abrahams transition rates in the form 



z/q exp 



2kT 



(5.3) 



with prefactor uq = Fi exp (— 26/q;), b being the distance between the neighboring 
sites on the conduction chain. Nearest-neighbor hopping within the model described 
by Eq. (|5.1|) and Eq. ( ^.3|) has been recently studied by computer simulations [|ble99 



Other forms of transition probabilities have also been considered in the literature. 
For example, a symmetric form was recently used in order to analyze analytically the 
field and temperature dependences of the hopping conductivity in one-dimensional 



systems ||dun96||: 



i±i,i = 1^0 exp 



2kT 



(5.4) 



Two models have been discussed in the literature with respect to the space- 
energy correlations of localized states responsible for transport. In the so-called 



Gaussian disorder model (GDM), such correlations are neglected |[bas93| , |ble99|| . In 



the correlated disorder model (CDM), it is assumed that energy distributions for spa- 
tially close sites are correlated, for example, due to dipole or quadrapole interactions 
||dun9q , |ga?95| , |nov98a| , |nov98b |. 

The drift mobility of carriers in organic disordered media usually demonstrates 
exponential dependences on temperature T and electric field E in the form 



/i oc exp 



To- ' 
T 



exp 




(5.5) 



where Tq and Eq are parameters. Much attention has been paid in recent years to the 
dependence of the mobility on electric field [|dun96| , |gar95| , |nov98a] , |nov98b|| and it has 
been claimed that in three-dimensional (3D) systems the observed field dependence 



5.2. ANALYTIC THEORY 



65 



in a wide range of electric fields can be explained only by the CDM. We check below 
whether this is true for a one-dimensional (ID) system. 

As far as the dependence of the mobility on temperature is concerned, the situa- 
tion is widely believed to have been cleared up many years ago. Computer simulations 
ch81 ] and analytical calculations | gru84 , mov86 | in 3D case were carried out giving 



the dependence 



/i oc exp 



(5.6) 



with C ~ 2/3. This expression is widely believed to be universal for the above model 
and it is often used to determine disorder parameter a of assumed Gaussian DOS in 
various materials from the experimental measurements ||3ch99a , Qem96 |. However, it 
is necessary to emphasize that there is no agreement among researchers concerning 
the magnitude of the coefficient C for one-dimensional systems. Ochse et al. ||och99a| | 
used the 3D value C = 2/3 for columnar systems with one-dimensional transport, 
while Bleyl et al. ||ble99|| obtained C ~ 0.9 in their computer simulations for ID case 



with asymmetric transition probabilities. The analytic calculations of Dunlap et al. 
dun96[| for symmetric probabilities predict C = 1 in ID case. 

In this chapter the temperature dependence in Eq. (|5.6|) will be studied in order 



to clarify this dependence for the GDM and the CDM. This chapter concentrates 
on the contribution of the author in that field, i.e. in computer simulations ||koh01 |. 
However, the work was done together with the analytical work performed by Dr. 
Harald Cordes ||cor01 



The analytical results will be summarized in section f>.'2\ and 
may occasionally be mentioned together with the computer simulation results as well. 
In section |5.3| the results of Monte Carlo simulations will be gathered. 



5.2 Analytic theory 

The steady-state drift velocity of carriers in one-dimensional periodic systems can 
be written exactly using the general solution found by Derrida ||der83|| . 



Nh 



N-l 
k=0 



fc,fc+i 



N-l 

n 



- fc+l,fc 



N-l i 



1+ E n 



=1 i=i 



- fc+j,fc+j-i 



(5.7) 



where N is the periodicity of the system and h is the distance between nearest neigh- 
bor sites. This formula was recently used by Dunlap et al. ||dun96|| to study the 
field dependence of the drift mobility in a one-dimensional system with symmetric 



transition rates [|dun96|| . The drift mobility is related to the steady-state velocity 



as 
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/i=-. (5.8) 

However in experiments, the drift mobility is usually obtained by the time-of-flight 
technique. Here charge carriers pass only once through a sample of finite thickness 
[ble99| , |och99a|| and drift mobility is calculated by dividing the sample length Nb by 



the mean transit time toTV and electric field E: 

Nb , , 

i^ = r^- (5-9) 

Murthy and Kehr ||mur90|| have derived the analytical expression for the mean transit 



time of carriers through a given chain of + 1 sites 

fc=o ^ fc=o ^ '=''=+1 i=k+\ j=k+i ^ i J+1 

This formula gives the transit time for a charge carrier that starts on site and 
finishes on site A^. The time is averaged over many transit times through the same 
chain with given values of transition probabilities Tij. 

One should be cautious with the apphcation of Eqs. ( |5.9| ) and ( [5. 101 ) at low electric 



fields. Even without an applied electric field carriers starting at site will pass through 
the system solely due to a diffusion process. At low fields, diffusion dominates the 
motion of carriers and it would be wrong to use Eq. (|5.9|). In such cases it would be 
more appropriate to estimate the mobility via the diffusion formula 

e b^N^ , , 

A* = 777;^—- 5.11 

which presumes the validity of the Einstein relation. 

Note that Eqs. ( ^.7|) , (|5.8| ) and Eqs. ( |5.9|) , (|5.10| ) are valid for any type of nearest- 



neighbor hopping of noninteracting carriers in one-dimensional systems. 

5.2.1 Drift mobility in the random-barrier model 

We start our analysis with a simple random-barrier model (REM). In the REM, all 
sites have equal energies being separated by energy barriers. Transition probabilities 
in the presence of electric field E are given by expressions 

/ ^^-\ebE^\^^-\ebE\ \ 
rfc,fc+i = I/O exp I ^ 1 (5.12) 



and 



= I/O exp I I, (5.13) 
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where is the energy barrier between sites k and k + 1. 

In the hmit of long periods (long chains), ^ 1, calculation of the drift mobility 
via Eqs. ( ^.7|) , (|5.8|) or Eqs. (|5.9| ), ( ^.101) for any chain with a given distribution 
of transition rates is equivalent to the averaging of Eq. ( |5.7| ) or Eq. ( |5.10| ) over all 
possible distributions of transition rates determined by disorder. In the following we 
will call this "averaging over disorder". Within the RBM, the averaging over disorder 
is given by ||der83| , |mur90 



k+l,k 
k,k+l I 



k,k+l 



(5.14) 



where (...) denotes the averaging over the distribution of barrier heights p(A). Note 
that Eq. (|5.14| ) can be derived by both equations, Eq. ( ^.7|) and Eq. ( ^.10|) . The 
explicit expressions for the averaging in the RBM are 
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one obtains for the drift mobility 
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(5.18) 

where erf(x) = ^ /q^ exp{—y'^)dy is the error function. In Fig. |5.1| the calculated field 
dependences of the drift mobility are shown for parameters b = 3.6 A, cr = 50 meV, kT 
= 25 meV. The solid line represents the exact solution for infinite chain given by Eq. 
( |5.18| ). All points in the figure correspond to mobilities for finite chains of = 500 
sites averaged over 1000 different chains. Circles were obtained via Eqs. ( |5.9| ), ( |5.10| ) 
while squares via Eqs. (|5.10| ), ( |5.11|) . Results of our Monte Carlo computer simulation 
are shown by crosses to demonstrate the excellent agreement of the simulation results 
with the analytical theory. At low fields, drift approximation (Eqs. ( p. 91 ), ( p.lOD ) 
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Figure 5.1: Field dependence of the carrier mobility in the RBM for Gaussian distribution 
of barrier heights with the width a = 50 meV at temperature kT = 25 meV. The solid line 
represents exact solution for the infinite chain. Data shown by circles and squares were 



calculated via Eq. ( 5.10 ) using drift and diffusion relation, respectively. Data shown by 
triangles were calculated via Eq. ( [5.7D , while those shown by crosses were obtained by 
Monte Carlo simulations presented in the section |5.3| . The chain length was 500 and 
averaging was performed over c = 1000 chains. 



for finite systems leads to the increase of the calculated mobility with decreasing 
field strength. Similar results are obtained for all considered models of disorder and 
various forms of the transition probabilities between neighboring sites. This result 
simply reflects the fact that charge carriers can penetrate through a finite system 
via diffusion motion even in the absence of an external electric field. By using Eq. 
( |5.9| ) one overestimates the mobility at low fields. The same happens with using Eq. 
( |5.11| ) at higher fields. Comparison of Eqs. ( |5.9| ) and ( p.ll| ) reveals the strength of the 
electric field at which a transition from the diffusion-controlled to the drift-controlled 
transit times takes place: E ^ 2kT/eNb. 

At low fields, transition time toAf does not depend on the strength of the elec- 
tric field E being determined mostly by diffusion. Using Eq. ( p.9| ) one artificially 
obtains an increasing drift mobility at decreasing field E. In numerous publications 
experimental results were reported that evidence the decrease of the drift mobility 
with increasing electric field ||pel88| , ^ch92| , |abk92|| . This was always observed at high 
temperatures at which diffusion can really dominate the motion of charge carriers. 
Experimental evidence for decreasing mobility with increasing field is usually obtained 
using equations similar to Eq. ( p.9| ), where the drift mobility is determined via the 
measured transit time by dividing the sample length by the value of this time and 
by the value of the field strength. We believe that this calculation is not appropri- 
ate and a diffusion equation similar to Eq. ( |5.11| ) should be used at low fields and 
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high temperatures. In 3D case, this equation should be shghtly modified, though the 
conclusion is the same: decreasing of the drift mobility with increasing field strength 
at low fields is an artifact caused by using the wrong equation [similar to Eq. ( |5.9|) ] 
for evaluation of the mobility on the basis of the measured transit time. In our ID 



calculations we illustrate this idea by using Eq. ( 5.11 ) instead of Eq. ( p.9| ) for finite 
systems at low electric fields. The result is shown by squares in Fig. |5.1| . The excellent 
agreement with the exact calculation for the infinite system at low fields confirms our 
conclusion. For the chosen parameters, the linear transport regime with the mobility 
independent of the electric field is valid up to the field strength of approximately 
10® V/cm, where a nonlinear regime starts with the mobility increasing with electric 
field. The latter regime is valid in a rather narrow range. At higher electric fields it 
is replaced by the trivial regime in which electric field eliminates the energy barriers 
between the sites of the chain and the transit time becomes field independent. In 
such a case mobility decreases proportionally to 1/E. 

In Fig. ^]2| the field dependence of the drift mobility is shown for a lower tem- 
perature kT = 15 meV for the same system as in Fig. |5.1| with a = 50 meV and b 
= 3.6 A. Clear mesoscopic effects are seen in this figure as evidenced by the differ- 
ence in the results obtained by averaging the mobilities (open symbols) and those 
obtained by averaging the transit times or the inverse mobilities (solid symbols). In 
all cases, averagings were performed over 1000 different chains, each chain consisting 
of 500 sites. The difference in the results of the two different averaging procedures 
refiects the presence of "fast" and "slow" chains. Averaging of times corresponds to 
the successive arrangement of chains where slow electron transitions over untypically 
high energy barriers determine the transport coefficients. On the contrary, the av- 
eraging of the mobilities corresponds to the parallel arrangement of chains, in which 
the fastest chains with untypically low energy barriers provide fast transit times of 
charge carriers. In the next section, we show that qualitatively the same effects are 
also inherent for the more realistic random-energy model. 



5.2.2 Drift mobility in the random-energy model without 
correlations 

In this section we consider a random-energy model with a Gaussian DOS described by 
Eq. (|5.1| ), presuming that there are no correlations between spatial positions of chain 
sites and their energies. This model is called in the literature a Gaussian disorder 
model (GDM). It has recently been suggested for the transport of charge carriers 
in discotic columnar liquid crystalline glasses with transition rates described by Eq. 
( |5.3| ) and it was studied by computer simulations ||ble99| ]. The exact analytic results 
for this model based on Eqs. (|5.71 ) and (|5.10|) will be presented below. 



In Eq. (p. 3D, site energies depend on the strength of the electric field E. They 
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Figure 5.2: Field dependence of the carrier mobility in the RBM. Chain length = 500, 
number of chains c = 1000, width of barrier height distribution a = 50 meV, temperature 
kT = 15 meV. Open circles and open triangles show the data obtained by averaging over 
chain mobilities calculated by Eq. ( 5.1C| ) and Eq. ( |5.7D , respectively. Solid symbols give the 



data obtained by the corresponding averaging over the inverse mobilities. 

are related to the zero field site energies 0^ as = ipk ~ ekbE. In accordance with 
Eq. (|5.3| ), the ratio of forward Tk,k+i and backward Tk+i,k hopping rates for any pair 
of neighboring sites is 

expf^^li^) (5.19) 



^k,k+l V kT 

equivalent to the condition of the detailed balance, which recently ||k:eh96|| has been 
used to determine the diffusion coefficient in thermal equilibrium for any form of the 
transition rates. 

Inserting relation ( |5.19| ) into Eq. (|5.7|), one obtains 
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(5.20) 
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fc=0 ^"•"^-^ k=0 ^ ' i=l exp[-^]^k+^,k+^+l 



The only random variables in Eq. (|5.20|) are the zero field site energies (pk and the 



forward transition rates Tk+i,k+i+i, which in turn depend only on neighboring site 
energies (pk+i and (pk+i+i- Since in all products appearing in the second term of the 
denominator, 0^ is not correlated with (p^+i or Fj^+j jt+j+i, averaging over disorder can 
be carried out for exp(— 0fe/A;T) and exp(0fc+j/A;T)F^^^ separately. Therefore 
the steady-state velocity and the mean transit time of carriers through long chains, 
3> 1 are given by 
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(5.21) 



(5.22) 



The relation v = bN/toN provides the equivalence of Eq. ( |5.22| ) and ( |5.21| ) in long 
chains as should be expected. Note that Eqs. ( |5.21| ) and ( |5.22| ) are exact for all 
kinds of nearest-neighbor hopping rates fulfilling Eq. ( |5.19|) . Thus they can as well 
be applied to the symmetric transition rates given by Eq. ( |5.4|) . 

The results of the averaging over Gaussian disorder for various terms in Eqs. 
(lOTI) , (lO^) are: 
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Thus in the one-dimensional GDM the drift mobility of carriers on an infinite chain 
is 
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At low electric fields this exact expression can be approximated by 
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while at high fields /i ~ v^b/ E. Interpolation of both approximations gives 



(5.27) 
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Field dependences of the drift mobility are shown in Fig. Q]3|for h = 3.6 A, cr = 50 meV, 
kT = 25 meV. The solid line represents the exact solution for an infinite chain given 
by Eq. ( p.26| ). Dashed lines show approximations via Eqs. ( |5.27| ) and fi ~ uob/E, 
respectively. The dotted line illustrates the interpolation formula [ Eq. ( ^.2^ )] which 
appears to be in excellent agreement with the exact solution. Circles and triangles in 
the figure show the calculated results for finite systems of 2000 sites averaged over 1000 
different chains once using the averaging of inverse transit times (circles) and once 
using the averaging of inverse velocities (triangles). Qualitatively, the results for the 



random-energy model in Fig. 5.3 resemble those for the random-barrier model from 



Figs. 5.1 and 5.2 and we refer the reader to the corresponding parts of the previous 



section, where the latter results were discussed. In Fig. p^ , the field dependences of 
the drift mobility are shown in the Poole- Frenkel representation (In /i versus ^/E) for 
two different temperatures and two different averaging procedures. The figure clearly 
shows that in the chosen model the analytic solution of the carrier mobility hardly can 
be described by the Poole- Frenkel law In /i oc ^/E. This conclusion is in agreement 
with previous studies of Gartstein and Conwell |^ar95|| , Dunlap et al. [ |dun96|| and 
of Novikov et al. |[nov98a| , |nov98b|| . Following these authors we consider in the next 
section the field dependence of the drift mobility in the model of correlated disorder 
(CDM). In the rest of this section we concentrate on the temperature dependence of 
the drift mobility at low electric fields in the Gaussian disorder model (GDM). 

The low field mobility for a finite chain is given by Eqs. ( |5.20 ) and (|5.8| ) in the 
hmit of E ^ 0: 
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For the infinite chain this expression is equal to 
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With hopping rates described by Eq. ( |5.3|) low field mobility in the GDM is exactly 
given by 
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Figure 5.3: Field dependence of the carrier mobility for nearest neighbor hopping with 
Miller-Abrahams rates between sites with Gaussian energy distribution. The solid line 
shows the exact solution for an infinite chain. Dashed lines correspond to the low-field 
and high-field approximations, while the dotted line gives the interpolation between them. 
Circles show the averaged mobilities of c = 1000 chains with = 2000 sites each calcu- 
lated according to Eq. ( 5.101 ) and Eq. ( |5.9| ). Triangles represent the results obtained by 
averaging the inverse mobilities for the same chains calculated by Eq. ( 5.2C| ) and Eq. (|5.8| ). 
Temperature and the width of the energy distribution were chosen as kT = 25 meV and a 
= 50 meV. 




Figure 5.4: Poole-Frenkel plots for carrier mobilities at kT = 25 meV (upper curves) 
and kT = 15 meV (lower curves). Circles show the averaged mobilities calculated by Eq. 
( |5.20D , squares show the corresponding results obtained by averaging of inverse mobilities. 
Number of chains was c = 10000 with N = 500 sites each. 
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Figure 5.5: Temperature dependence of the low-field mobility for cr = 50 meV. Solid line 
represents the solution for the infinite chain given in Eq. ( 5.33| ). Circles and squares show 
the results obtained by averaging of mobilities and averaging of inverse mobilities calculated 
by Eq. ( ^.31 ), respectively, with averaging over c = 1000 chains of iV = 2000 sites. Upward 
and downward triangles are the corresponding values for c = 10000 and A'^ = 200. The 
dashed line indicates for illustration the temperature dependence with the slope equal to 
-1. 



This equation is valid for the infinite chain. The corresponding temperature depen- 
dence of the low-field mobility is shown by the solid line in Fig. |5.5| . The slope of 
Infi vs. {a/kTy in this curve differs slightly from —1, due to the temperature de- 
pendence of the preexponential factor in Eq. ( ^.31| ). This result is in good agreement 
with the recent computer simulations of Bleyl et al. [|ble99|| and with the analytic 
calculations of Dunlap et al. ||dun96|| , although the latter analytic calculations have 



been carried out for correlated systems with symmetric transition rates described by 
Eq. ( |5.4| ). This shows that in ID systems Eq. ( ^.6| ) with C ~ 1 is correct and even 
stable against the choice of the form of transition rates. 

The temperature dependences of the low-field drift mobility in finite systems 
calculated according to Eq. ( 5.29| ) are also shown in Fig. |5.5| . These results are 
really striking. They clearly show that the temperature dependence of the mobility is 
related to the size of the system. This is a manifestation of the mesoscopic character 
of the hopping transport problem which is mostly pronounced in ID systems. As 
has been already discussed with respect to the RBM, the averaging of transit times 
or inverse mobilities over various finite chains corresponds to the calculation of the 
mobility in a long system consisting of all these chains connected sequentially one after 
another. On the contrary, the averaging of the mobilities over various finite chains 
corresponds to the calculation of the carrier mobility in the system in which these 
chains are arranged in parallel to each other. In the latter case, the "fast" chains with 
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untj^ically short transit times strongly contribute to the averaged value of the carrier 
mobility and they dominate the transport of carriers. This is the very essence of the 
mesoscopic effects [|rai91j] . Having in mind an application of our results to discotic 
organic disordered systems where many quasi-lD current channels are connected in 
parallel being perpendicular to the sample surface, one should conclude that the 
temperature dependence of the electrical current in such systems does change with 
the thickness of the samples. For example, a comparison between the data obtained 
by the averaging of mobility values for chains with = 2000 sites and chains with 
= 200 sites suggests that for shorter chains and hence for the thinner samples 
the temperature dependence of the low-field drift mobility should be weaker than 
that for thicker samples. Of course, it would be desirable to verify our prediction 
exp er iment ally. 

We now turn to the question how decisive for our results is the choice of the 
Miller- Abrahams transition probabilities described by Eq. (|5.3| ) and used for the above 
calculations. For comparison, one can insert into described mathematical apparatus 
symmetrical transition probabilities described by Eq. ( p.4| ). One straightforwardly 
comes via Eqs. ( ^.21|) and ( |5.8| ) to the result 
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(5.32) 



In Fig. we compare the field dependence of the mobility given by Eq. ( |5.32| ) 
with that for Miller-Abrahams rates. After an appropriate scaling of the results 
for symmetric rates the mobilities in both cases (symmetric and Miller-Abrahams 
hopping probabilities) are the same at low electric fields. This confirms the conclusion 
of Dunlap et al. ||dun96|| that the form of transition probabilities is not essential for 
the field dependence of the drift mobility at low fields, provided the detailed balance 
is fulfilled. Dunlap et al. [|dun96|| obtained this result for correlated systems (CDM), 
while we confirm their conclusion for the uncorrelated Gaussian model (GDM). In 
the non-linear regime at higher fields {E ~ 10^ V/cm) deviations between the results 
for symmetric and for Miller- Abrahams hopping probabilities are seen in Fig. B^. 



5.2.3 Drift mobility in the random-energy model with corre- 
lated disorder 

So far we have considered the hopping transport of charge carriers along regular chains 
of sites with Gaussian distribution of site energies and without correlations between 
the spatial positions of sites and their energies. Such a model has recently been 
suggested by Bleyl et al. |[ble99|| for transport of charge carriers in columnar discotic 



liquid-crystalline glasses. Bleyl et al. treated this model by computer simulations 
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Figure 5.6: Influence of hopping rates on tlie field dependence of the carrier mobility. 



Circles show the results calculated by Eq. (5.20) with Miller- Abrahams rates defined in 



Eq. (5.3). Squares display the mobilities calculated with symmetric hopping rates defined 
in Eq. (5^). Above data were obtained for 1000 chains each containing 500 sites. Solid 
lines show the corresponding analytic solutions for infinite chains. For a better comparison, 
triangles are plotted indicating the mobilities with symmetric hopping rates scaled by factor 
iexp[cTV(2A:r)2]. 



and obtained results similar to our exact analytic results from the previous section. 
We then analyzed whether or not the assumption on the uncorrelated distribution of 
site energies is essential for these results. Such a question was raised by Gartstein and 
Conwell |^ar95|| , by Dunlap et al. ||dun96|| and by Novikov et al. ||nov98a| , |nov98b 



Using a model similar to that of Gartstein and Conwell, we present below, exact 
analytic solutions for the drift mobility of charge carriers in systems with correlated 
disorder (CDM). 

In order to construct a model with correlated disorder (CDM) we first generate 
provisional site energies ipk- The zero-field energy 0^ of a charge carrier on site k is 
then determined by the arithmetic average of provisional energies ip of neighboring 
sites: 



0fc = ^ ^ H '^k+j (5.33) 

j=-m 

Here A = 2m + 1 is the correlation length in units of the intersite separation h. To 
ensure that the resulting site energies 0^ have a Gaussian distribution with variance 
a we distribute energies ifj according to 
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Due to the correlations, two site energies at zero field (f)k and (f)k+i are not independent 
as long as i < A. This prevents the possibility to carry out the averaging over disorder 
in the style as it was done above on the way from Eq. ( ^.20|) to Eq. (|5.21|) . Instead 
we have to split the sum in the second term of the denominator in Eq. ( [5. 201 ) into 
two terms: 
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(5.35) 
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In the second term of the denominator the site energies (pk are not correlated with 
site energies (pk+i or with transition rates T^+i^k+i+i- Hence here the averaging of 
exp(— 0fc//cT) and [exp{—(j)k+i/kT)Tk+i,k+i+i]~^ over disorder can be carried out sep- 
arately. Inserting Eq. ( p.33| ) into the last term of the denominator one obtains 
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With our choice of site energy correlation, Eq. ( |5.33| ), the energy difference between 
neighboring sites k + i and k + i + 1 is {ipk+i+i+m — i^k+i-m)- The transition rates 
Tk+i,k+i+i are uncorrelated with ipk+m+j and ipk-m+i+j for all j and the averaging for 
the corresponding terms can be performed independently. Note that for long-range 
correlations one has to average the entire product in Eq. ( ^.36| ). In the latter case 
Eq. ( |5.35| ) cannot be used and one has to solve Eq. ( [5. 201 ) instead. This may prevent 
an exact solution. For correlations given by Eq. ( |5.33| ) we get 
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Averaging over disorder with using Gaussian distribution of site energies g{ip) given 
by Eq. ( ^.34|) leads to the following result for the drift mobility of carriers in an infinite 
chain within the CDM 
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In Fig. |5.7| the corresponding field dependences are plotted for different correlation 
lengths A for finite chains (symbols) as calculated via Eqs. (|5.2CI|) , ( pl8|) and for the 
infinite chain [Eq. ( |5.38|) ]. The particular shape of the field dependence of the drift 
mobility in the CDM for the infinite system has been already discussed in detail 
by Gartstein and Conwell |^ar95|| , by Dunlap et al. [|dun96|| and by Novikov et al. 
|nov984 |nov98b|| . We would like to focus our attention here on the other aspect of 



the phenomenon, namely, on the increasing difference between the results obtained 
by averaging of mobilities and those obtained by averaging of inverse mobilities in 
finite systems with increasing correlation length. This trend is clearly related to the 
smaller number of high barriers for charge carriers in systems with longer correlation 
lengths. The mesoscopic effects become more pronounced in systems with longer 
correlation length. The increase of the correlation length in the CDM is analogous to 
the decrease of the total number of sites in the chain in the GDM. 



5.2.4 Concluding remarks 

Analytic calculations were presented for the hopping drift mobility of charge carriers 
in a one-dimensional regular chain of localization sites with energetic disorder. Exact 
results were obtained for temperature and field dependences of the mobility in systems 
with uncorrelated disorder (GDM) and those with a correlated disorder (CDM) for 
both asymmetrical and symmetrical transition probabilities. These dependences are 
shown to be determined by the length of the system. Therefore even the exact 
solutions usually obtained for infinite systems might be of low value for thin samples 
which are used in real devices. Focusing on the description of the transport properties 
of the discotic liquid-crystalline glasses, one should take into account that the real 
device samples usually contain only 50 to 100 sites ||chr97b||. For such systems one 



should use the theoretical results obtained in this report for finite systems and not 
those derived for infinite chains. 

In this section we considered only transitions of charge carriers to the nearest 
sites on a conducting chain. In the next section (section p.3| ) it is shown by a Monte 
Carlo computer simulation that in a correlated system (CDM) transitions to more 
distant neighbors do not play an essential role while for uncorrelated systems such 
distant transitions can be important. Hence in the CDM our analytical results can 
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Figure 5.7: Influence of the space-energy correlations on the field dependence of the 
carrier mobility at kT = 15 meV and a = 50 meV. The different correlation lengths are 
A = 1 (circles), A = 3 (squares), A = 11 (upward triangles) and A = 101 (downward 
triangles). Open symbols show the averaged mobilities calculated by Eq. ( |5.20| ). Solid 
symbols show the corresponding results obtained by averaging of the inverse mobilities for 
(A^ = 1000, c = 1000). Solid lines show the solutions for the infinite chain given by Eq. 



be definitely considered as exact. In the case of uncorrelated disorder (GDM), longer 
electron hops can modify the results for the field and the temperature dependences 
of the carrier drift mobility as shown in section |0. 



5.3 Monte Carlo simulations 



In section ^]2| analytic calculations were carried out for the steady-state drift mobil- 
ity of charge carriers in one- dimensional systems with taking into account hopping 
transitions only between the nearest localized states. Special attention was paid to 
the temperature and field dependences of the mobility for various hopping models. 
We check here the corresponding results by a Monte Carlo computer simulation and 
show that analytic theory works perfectly. However, this theory is unable to take into 
account hopping events to localized states situated further in space than the nearest- 
neighboring sites. We will show below that electron transitions to more distant sites 
decisively influence the most fundamental transport properties, such as, for example, 
the temperature dependence of the drift mobility. This dependence at low electric 
fields is described by Eq. (^.6|). Such a dependence is observed experimentally and 
it has also been obtained in computer simulations and analytic calculations. In the 
three-dimensional (3D) case, coefficient C in Eq. ( |5.6|) is believed to have the magni- 
tude ||bas93 1 C ~ 2/3, while in the one-dimensional (ID) case the value C ~ 1 has 
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been obtained for symmetric transition rates determined by Eq. (|5.4| ) ||dun96|| . Equa- 
tion (|5.6|) is widely used to determine parameter a of the density of states (DOS) in Eq. 
(|5.1|) for various materials from experimental measurements of the ln(/i) vs. (1/T)^ 
dependences, in particular for the ID hopping transport | |och99a| , |ble99| | . Therefore, 
the magnitude of the coefficient C is really important. This coefficient C in Eq. (|5.6|) 
cannot be universal. Equation ( p.2|) shows that both the energy difference between 
localized states involved in the carrier jump and the distance between these states 
determine the hopping probability. In such processes, characteristic temperature Tq 
should depend on the concentration of localized states N and on the decay length 
of the carrier wave function in the localized states a in the combination Na in the 
ID case and Na^ in the 3D case ||shk84| ] . In the computer simulations ||sch81|| this 



parameter was taken as Na^ = 0.001 and it is not surprising that some particular 
value for coefficient C was obtained. A very close value C ~ 2/3 was also obtained 
by the effective-medium theory in the 3D case ||gru84| , |mov86|| . It is not correct to use 
this value for ID systems. Concentrating on the ID transport we emphasize that the 
value of the coefficient C = 1 obtained in the analytic calculation for Miller- Abrahams 
transition rates and also known from the literature for symmetric rates ||dun96|| can 
also be wrong. This value was obtained only for the nearest-neighbor hopping. Be- 
low we show that in systems with uncorrelated energetic and spatial distributions of 
localized states, hopping transitions to more distant sites than the nearest neighbors 
influence the value of the coefficient C, which is decisive for determination of the dis- 
order parameter a of the DOS for various materials from experimental measurements 
of the ln(/i) vs. (l/T)^. 



5.3.1 Numerical algorithm and simulation details 



In our simulations a Gaussian density of states (DOS) was assumed following the data 
from the literature ||bas93| , |dun96| , |och994 |ble99|| . In the so-called Gaussian disorder 
model (GDM), the energy distribution was described by Eq. (|5.1|) and no correlations 
between energy and space positions of localized states were allowed. In the so-called 
correlated disorder model (GDM) ||dun96| , |gar95| , |nov98a|| , the site energies were cor- 
related with energies of neighboring sites. We used correlations described by Eqs. 



( |5.33| , |5.34|) with m equal to 1. In the simulation procedure, c independent regular 
chains of localization sites, each chain containing sites were studied simultaneously. 
Only one particle was considered on each chain. The simulation algorithm was the 
following. A charge carrier starts its motion on the beginning site of a chain (site 
number 1). The simulations were stopped after 90% of carriers reached the final sites 
of the chains (sites number A^). The rate Fjj of a hopping transition from an occupied 
site i to an empty site j were determined by Miller- Abrahams equation (Eq. |5.2|) . The 
lifetime of a carrier at a site i was calculated as 



5.3. MONTE CARLO SIMULATIONS 



81 



t. = ^ (5,39) 

2^ J- ij 

j 

where x is a random number from the uniform random-number distribution between 
and 1. The sum in Eq. ( [5.39| ) was taken over sites with numbers j = i ± 1 for 
nearest-neighbor hopping and j = i ± 1, i ± 2, i ± 3, i ± 4, i ± 5, i ± 6 for simulations 
with tunneling to further neighbors. 

For the fastest 10% of particles two different methods to calculate the mobility 
were used. In the first method, transition times for different chains were averaged 
and the obtained value (t)io% was used to calculate the carrier mobility via relation 
/^io% ~ E(t)u-,<y ' '^^^^^ d is the length of the chains and E is the strength of the applied 
electric field. The other method was based on the calculation of carrier mobilities for 
the fastest 10% of carriers and on averaging the mobility values. The latter method 
in our simulations is equivalent to the averaging of the inverse transit times (l/t)io%. 
The averaged mobility in this method was calculated as = . The same 

averaging procedures were performed for the 90% of charge carriers with finding the 
corresponding quantities = e^^^, and /x^o?,f^ = 

Simulations were carried out for various magnitudes of the electric field E and 
temperature T. In order to obtain the value of the zero-field mobility fizf, the pro- 
cedure suggested by Bleyl et al. ||ble99|| has been used. The data were plotted in the 



Poole-Frenkel (PF) representation 

ln(/i) = ln(/i,/) + (5.40) 
with the aim of finding via extrapolation to zero field. To illustrate this extrap- 



olation procedure, we show some of our results in Fig. |5.8| . These particular results 
were obtained for the GDM with the energy scale of the DOS a = 30 meV and the 
temperature range corresponds to kT between 10 meV and 30 meV. The mobility val- 
ues were obtained by averaging of the inverse transit times for 90% of charge carriers. 
After having found /i^/ for various temperatures, the T dependence of this quantity 
was plotted as 



/^^/ l^j =/^oexp 



(5.41) 



where /xq is a temperature-independent prefactor and C is a constant to be determined 
from such plots. 

Obtained temperature dependence described by Eq. ( |5.41| ) is plotted in Fig. ^ 



illustrating some of our data obtained for the nearest-neighbor hopping in the GDM. 
A similar Monte Carlo simulation was carried out recently by Bleyl et al. ||ble99|| . 



who obtained very similar results. Our data confirm the coefficient C ~ 0.9 in Eq. 
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Figure 5.8: Field dependences of the carrier mobility at various temperatures in the 
Poole-Frenkel representation for the GDM with o" = 30 meV. 



(|5.41|) obtained earlier by Bleyl et al. However, computer simulations of the nearest- 
neighbor hopping are not of considerable interest, as exact results for such processes 
can be obtained analytically for finite and infinite systems as previously shown in 
the analytical work (section |5.2D . Hopping processes to further neighbors than the 
nearest ones can hardly be treated exactly by analytic theories. In this section we use 
computer simulations to study whether such distant hopping transitions are important 
for transport coefficients, like the drift mobility, or not. 

In order to answer this question, we performed computer simulations of hopping 
processes to distant neighbors for both GDM and GDM and compared the results with 
those obtained only for hopping to the nearest neighbors. The results for the latter are 
given in section |5.3.2| , while the results for distant hops are presented in section |5.3.3 



Also the question will be addressed on the significance of the space-energy correlations 
for transport processes as given by the comparison of the results for the GDM and 
GDM. The latter question has recently been studied with respect to the dependence 
of the carrier drift mobility on the electric field ||dun96| , |gar95| , |nov98a| , [nov98b|| and 
it was shown that the space-energy correlations play an essential role for the field 
dependence of the mobility. Below we mainly consider the temperature dependence 
of the carrier mobility. The quantity of our main interest will be the coefficient C in 
Eq. ( [5.41| ). Before presenting the results, we would like to note that the extrapolation 
procedure according to Eq. ( p.40| ) cannot be exact. We have found that the results 
for fi^f can slightly depend (by less than 10%) on the chosen range of field strengths 
used for extrapolation. We always used the interval 9 x 10^ < -E < 2.5 x 10^ V/cm 
where the PF dependence seems to be valid as shown in Fig. |5.8|. Moreover, the values 
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Figure 5.9: Temperature dependence of the zero-field drift mobility according to Eq. 
( |5.41 ) for the GDM with o" = 30 meV (diamonds) and o" = 50 meV (crosses). Dashed and 
dotted lines represent Eq. ( 5.41| ) with C = 0.928 and C = 0.924 respectively. 



of the determined coefficient C (by less than 10%) depend shghtly on the way how 
parameter a/kT was changed - by changing T with fixed a, or by changing a at fixed 
T. This is caused by the extrapolation procedure based on Eq. ( |5.40|) , because T and 
a determine the field dependence of the carrier drift mobility not in the combination 
a/kT which is essential for the temperature dependence of the zero- field mobility as 
has been clearly shown by Dunlap et al. ||dun96|| . Therefore, not the absolute values 
of the coefficient C, but its trends with changing the simulation parameters will be 
of most interest with respect to our simulation results. 



5.3.2 Nearest neighbor hopping 

All simulation results given below were averaged over c = 1000 chains, each containing 
= 2000 sites. Following Bleyl et al. [ ble99| ] only Miller- Abrahams transition rates 
determined by Eq. ( |5.3| ) were used. Parameters Fi = 1.875 ■ 10^^s~\ b = 3.6 A, a 
= 1 A were chosen, where b is the distance between the neighboring sites. As long 
as only the nearest-neighbor hopping is considered, the choice of parameters b and 
a cannot infiuence the temperature dependence of the mobility which is of prime 
interest in our study. In order to illustrate the simulation results, we plot in Fig. 
5.10| the field dependence of the drift mobility at kT = 25 meV and cr = 50 meV for 
both GDM and CDM. Results of our computer simulation are shown by crosses for 
the CDM and by diamonds for the GDM. Dashed lines represent the results of the 
exact analytic calculations for the particular chains used in the simulation procedure 
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based on the methods described in section |5.2| . Excellent agreement with the exact 
analytic results confirms the validity of our simulation procedure. In agreement with 
previous calculations ||dun96| , |gar95|, |nov98a|] , one can see in Fig. |5.1CI| that the mobility 



values are higher and the field dependence starts at lower fields for the CDM than 
for the GDM. We also present in Fig. 5.1CI| the statistics for the average number of 



hops necessary for a carrier to penetrate through a chain containing 2000 sites. At 
low fields, the motion is diffusive and the number of hops is enormously larger than 
the number of sites indicating that charge carriers often move back and forth. This 
supports the conclusion of analytical calculations (section p.2\) that at low fields one 
should use the diffusion equation [Eq. ( |5.11| )] rather than the drift approach [Eq. 



(|5.9|)]. Erroneous application of Eq. ( p.9| ) at low fields leads to the artificial increase 
of mobility values with decreasing field strength. Unfortunately, this artifact is often 
treated as a true physical effect. At very high fields, the carriers mostly hop in the 
direction prescribed by the field as indicated by the equality between the number 
of hops necessary for charge carrier to penetrate through the chain and the number 
of chain sites (A^ = 2000). In such extreme conditions, the drift velocity does not 
depend on the field strength and Eq. ( ^.8|) leads to the relation n oc 1/E for high 
fields. Performing the extrapolation and fitting procedure described in section |5.3.1| , 
we obtained the values for the coefficient C in various conditions. These values 



were obtained by the four different averaging methods described in section |5.3.1 
In the following, CJo% example, will denote the coefficient C corresponding to 
fastest 10% of charge carriers (i.e., "fastest" 10% chains) and "tav" means that the 
values were obtained by averaging of arrival times. Values for the coefficients 
were correspondingly obtained via the averaging of mobilities. Obtained results are 
gathered in Table ^.1| , where different horizontal lines correspond to simulation results 
at different values of the disorder parameter a in Eq. ( |5.1|) . Temperature values were 
changed so that parameter a/kT had magnitudes between 1 and 3. In all cases, 
the values of the coefficient C are between 0.9 and 1, in good agreement with the 
results of previous computer simulations |[ble99|| and analytical calculations ||dun9(: 



(see also the analytical solution in section |5.2|) . No striking difference concerning 
the temperature dependence of the carrier drift mobility has been found between the 
GDM and CDM as can clearly be seen from Table f).l\ 

5.3.3 Hopping with tunneling to further neighbors 

In the previous section as well as in the analytical calculations (section |5.2|) , only 
transitions of charge carriers to the nearest-neighboring sites were considered. In 
such conditions, transition rates are determined by Eq. ( |5.3| ). Now we would like to 
check whether and at what conditions, if at all, the tunneling hops to further sites 
than the nearest neighbors influence the transport coefficients. In order to answer 
this question we repeated all simulations described in section p.3.2|, but this time al- 
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Figure 5.10: Left axis: dependences of the carrier mobility ^qq^ on the electric field for 
the GDM and CDM. The lines show the results of the exact analytic solutions for system 
parameters a = 50 meV and kT = 25 meV; right axis: average numbers of forward steps as 
functions of electric field. The chain length was equal to 2000 sites in all simulation runs. 



Table 5.1: Coefficients C in Eq. (|5.41| ) for the nearest-neighbor hopping 



model 


parameter 


•-"loro 


r^mobav 


r^tav 
'-"90% 


r^mobav 
'-"90% 


GDM 


a = 30meV 


0.89 


0.89 


0.95 


0.93 


a = AOmeV 


0.90 


0.89 


0.95 


0.93 


a = bOmeV 


0.90 


0.90 


0.94 


0.92 


CDM 


a = 28.867meV 


0.92 


0.92 


0.99 


0.96 


a = AOmeV 


0.92 


0.92 


0.98 


0.96 


a = 50meV 


0.91 


0.91 


0.96 


0.94 
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lowing the tunneling of charge carriers to further sites. For such distant hops explicit 
i?— dependent transition probabilities described by Eq. ( ^.2|) should be taken into 
account instead of the simplified probabilities described by Eq. (|5.3|) . Moreover, we 
carried out all simulations with changing the values of parameter a which determines 
the decay length of the carrier wave function in the localized states. The intersite 
distance was kept constant b = 3.6 A, the values of a were changed between 1 A 
and 3 A. Transitions to six neighbors in each direction were allowed in the simulation 
procedure. We have checked that tunneling to further states than these 12 neighbors 
does not play any role at all sets of parameters used in the simulation runs. Results 
obtained for coefficient C in Eq. ( |5.41| ) at various values of a and a and for various av- 



eraging procedures are gathered in Table 5.2. The main conclusion one can draw from 



this table is the following. For the CDM, i.e., for a system with correlated disorder, 
tunneling to further sites than the nearest neighbors do not play any essential role, 
while for the GDM, i.e., for uncorrelated systems these distant hopping transitions 
slightly influence the transport coefficients. The insignificance of distant transitions 
for systems with correlated disorder is unsurprising. In such systems, energies of 
neighboring sites are close to each other due to correlation effects. Thus the variable- 
range hopping is not signiflcant for such systems, because the closest in energy and 
space states are chosen by carriers in their motion and these states are the nearest 
neighbors. The quantitative conflrmation of this qualitatively clear statement is of 
great importance because it implies that for systems with correlated disorder one can 
use analytic theories described in section ^]2| and one can be sure that the restriction 
of the analytic theory which considers only the nearest-neighbor hopping is not at all 
essential for the obtained results. 

For systems with uncorrelated disorder the situation is also rather optimistic with 
respect to the exact analytic solutions (see section |5.2|) . From Table f).2\ it can be 



seen that even for rather high value of the decay constant a = 2 A which makes the 
distant hops very favorable, the value of coefficient C ~ 0.8 is not that dissimilar 
from its value C ~ 0.9 for the nearest-neighbor hopping. Thus even for systems with 



uncorrelated disorder, exact analytic solutions (see section |5^ ) taking into account 



only the nearest-neighbor hopping give reasonable values for transport coefficients. 
5.3.4 Concluding remarks 

Computer simulations of charge carrier transport in disordered ID systems show that 
for systems with correlated space-energy disorder the variable-range hopping does not 
play an essential role, implying that exact analytic solutions (section |5.2| ) based on 
the nearest-neighbor hopping are valid for such systems without any restrictions. For 
systems with uncorrelated disorder the analytic results for transport coefficients from 



section 5.2 can be considered as an approximation. For example, taking into account 



the more distant hops than hops to the nearest neighbors one obtains the value of 
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Table 5.2: Coefficients C in Eq. (|5.41|) when taking into account further hopping 
transitions than those to the nearest-neighboring sites. 



model 


parameter 


«(A) 


•-"10% 


/^mobav 
•-"10% 


r^tav 
•-"90% 


/^mobav 
•-"90% 






1.0 


0.86 


0.86 


0.92 


0.89 




a = SOmeV 


2.0 


0.79 


0.79 


0.86 


0.82 






3.0 


0.75 


0.75 


0.83 


0.79 






1.0 


0.87 


0.87 


0.92 


0.90 


GDM 


a = AOmeV 


2.0 


0.79 


0.79 


0.84 


0.82 






3.0 


0.74 


0.74 


0.79 


0.77 






1.0 


0.85 


0.84 


0.91 


0.87 




a = 50meV 


2.0 


0.77 


0.77 


0.82 


0.80 






3.0 


0.73 


0.73 


0.78 


0.76 






1.0 


0.93 


0.93 


1.00 


0.97 




a = 28.867me\/ 


2.0 


0.91 


0.91 


0.98 


0.96 






3.0 


0.88 


0.88 


0.96 


0.93 






1.0 


0.92 


0.92 


0.98 


0.95 


CDM 


a = AOmeV 


2.0 


0.90 


0.89 


0.96 


0.94 






3.0 


0.87 


0.87 


0.93 


0.91 






1.0 


0.92 


0.92 


0.97 


0.95 




a = bOmeV 


2.0 


0.89 


0.89 


0.95 


0.93 






3.0 


0.87 


0.87 


0.92 


0.90 
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the coefficient C in Eq. (|5.41| ) by approximately 10% different to that for only the 
nearest-neighbor hopping. With such an accuracy one can even use analytic results 
from section |5.2| for systems with uncorrelated disorder. 
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Summary 



The study of different physical properties of disordered systems by means of computer 
simulations was the aim of the present thesis. The incredibly huge and fast growing 
advances in information technology in the last decades, has allowed physicists to use 
computer resources more effectively. The extensive use of computers has enabled 
the development of a new and extensive part of physics, lying somewhere between 
experimental and theoretical physics. 

In order to describe the physical properties of amorphous materials it is very 
important that one has knowledge of the atomic arrangement. The structure of dis- 
ordered solids usually depends on their preparation method. It is therefore important 
that the preparation method is modeled first, following this, the received structures 
need to be studied. In chapter]^, low energy molecular dynamics simulations of atomic 
beam growth of amorphous carbon on a diamond [111] surface were carried out. Our 
aim was to construct large scale amorphous carbon models over 100 atoms by the 
quantum-mechanical treatment of the interatomic potential. The interactions were 
described by a well-tested and widely applied tight-binding potential. The structures 
were grown by atom-by-atom deposition onto substrates. The growth process was 
simulated as in experiments without any artificial model of energy dissipation. We 
used two different average bombarding energies Ef,eam = 1 and 5 eV, and two differ- 
ent substrate temperatures T^^j^ = 100 and 300 K for 30 ps simulation time. The 
simulations were prolonged by 15 ps at Tsub = 100 K to obtain larger structures. Six 
networks were prepared by the deposition technique with periodic boundary condi- 
tions in two dimensions. The most important structural parameters, like densities, 
radial distribution functions, coordination numbers, bond angle distributions, and 
ring statistics were then analyzed. The deposition rate was much higher than it is 
in experiments. On the basis of our time development investigation we are certain 
that much lower deposit rates would not cause a significant difference in our mod- 
els. Our lower substrate temperatures, compared to experiments, provide a quicker 
energy dissipation which slightly compensates for the high deposition rate. 

In chapter ^, we simulated the atomistic simulation of the bombardment process 
during the bias enhanced nucleation (BEN) phase of diamond chemical vapor deposi- 
tion (CVD). Special attention was paid to the study of different hydrocarbon (C^^IIy) 
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projectile bombardment at varying impact energies. The simulations showed that 
larger kinetic energies produce deeper subplantations and that it is more probable 
that the species will be dissociated. For acetylene (C2H2) the carbon-carbon bond is 
almost never broken at low energies (20, 40 eV) and it is broken in 50 to 90% of cases 
with larger energies (80 eV). The computer experiments manifested that approxi- 
mately 60% of the projectile hydrogen-carbon bonds are broken during the subplan- 
tation. Usually, one of the hydrogen-carbon bonds is not broken for methyl (CH3) and 
the other two hydrogens form H2 molecules. At the typical 200 V bias voltage, which 
is optimal for diamond nuclei formation in BEN, the C2H2 projectiles have 80 eV, 
and CH3 ion species have 40 eV kinetic energy. There are significant differences in the 
structural rearrangements in the film caused by these two projectiles at this voltage. 
For the C2H2 species, the penetration is deeper and these projectiles cause changes 
in the deeper region than is the case with the CH3 species. The CH3 projectiles do 
not contribute to the deep variations in the film and cause structural rearrangements 
only in the surface regions. The structural changes in the film showed that higher 
kinetic energy ionic species under the average deposition depth generate a thicker 
sp^ rich film and the surface mostly consists of atoms with sp^ hybrids. The border 
between these two regimes move downwards with larger deposition energies. The 
results support the subplantation model. In the future it would be a very interesting 
task to do simulations over a silicon substrate and study the effect of this crystalline 
layer on the formation of diamond nuclei in the amorphous carbon film. 

Chapter ^ deals with the simulation of amorphous silicon structures. The same 
atomic deposition program code, which had been developed for studying the amor- 
phous carbon structure growth and was reported in chapter 0, was applied here for 
silicon. We used a well tested tight-binding Hamiltonian to describe the interatomic 
interaction between silicon atoms. A similar tight-binding potential was developed by 
the same group for carbon systems and we successfully applied it for the description 
of the amorphous carbon growth (see chapter |^). Nine different amorphous silicon 
structures were obtained by our deposition computer method. An additional tenth 
model was prepared by rapid cooling in order to check the difference between atom- 
by-atom deposition on substrate and melt- quenching preparation techniques. The 
most significant difference we have found between the models prepared in different 
conditions, is in the ring statistics, i.e. in the medium range order. Triangles are also 
present in the atomic arrangement, supporting three other independent pieces of evi- 
dence into the existence of triangles in amorphous silicon structures. Nevertheless, the 
fragments like equilateral triangles have never been considered at electronic density of 
states calculations or at an unsolved problems such as the breaking of weak bonds by 
prolonged illumination, which is the major mechanism for the light-induced defect cre- 
ation. The tight-binding model potential used in our computer experiments slightly 
overestimates the bond distances, by giving values which were observed previously in 
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liquid silicon forms. The potential gives overcoordinated atomic sites (fivefold, and 
even sixfold coordinated atoms) during the amorphous silicon structure simulations. 
The sixfold coordinated atoms can be found in the liquid phase of silicon instead of in 
the amorphous structure. The well-known Stillinger- Weber classical empirical poten- 
tial for amorphous silicon had a similar difficulty and has recently been successfully 
modified. A possible future aim would be to improve this tight-binding potential by 
a similar process. 

Finally, in chapter ^, calculations for the hopping drift mobility of charge carriers 
in a one-dimensional regular chain of localization sites with energetic disorder were 
performed. In the analytic calculations, exact results were obtained for temperature 
and field dependences of the mobility in systems with uncorrected disorder (Gaussian 
Disordered Model) and those with a correlated disorder (Correlated Disorder Model) 
for both asymmetrical and symmetrical transition probabilities. These dependences 
are shown to be determined by the length of the system. Therefore even the exact 
solutions usually obtained for infinite systems might be insignificant for thin samples 
which are used in real devices. Focusing on the description of the transport properties 
of the discotic liquid-crystalline glasses, one should take into account that the real 
device samples usually contain only 50 to 100 sites. For such systems one should 
use the theoretical results obtained in section |5.2| for finite systems and not those 
derived for infinite chains. In analytic calculations we considered only transitions 
of charge carriers to the nearest sites on a conducting chain. The results of Monte 
Carlo computer simulation reveal that for systems with uncorrelated space-energy 
disorder, the variable-range hopping plays an important role. Longer electron hops 
can modify the results for the field and temperature dependences of the carrier drift 
mobility as presented in section |573| . For such systems the analytic results for transport 
coefficients can be considered as an approximation. Taking into account, for example, 
the more distant hops than hops to the nearest neighbors one obtains the value of the 
coefficient C, which describes the temperature dependence of mobility value in Eq. 
( |5.41|) , which is approximately 10% different to that for only the nearest-neighbor 
hopping. With such accuracy one can even use analytic results (from section |5.2|) 
for systems with an uncorrelated disorder. However, in the case of correlated space- 
energy disorder, this effect is much lower, implying that exact analytic solutions based 
on the nearest-neighbor hopping are valid for such systems without any restrictions. 
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Ziel dieser Dissertation war die Untersuchung verschiedener physikalischer Eigen- 
schaften von ungeordneten Systemen mit Hilfe von Computer Simulationen. Die 
iiberwaltigend schnell gewachsenen Entwicklungen in der Informationstechnologie in 
den letzen Jahrzehnten hat es ermoglicht, Computer Ressourcen effektiver als je zuvor 
zu nutzen. Die ausgdehnte Nutzung von Computern hat eine neues und umfangre- 
iches Gebiet der Physik hervorgebracht, das irgendwo zwischen experimenteller und 
theoretischer Physik hegt. 

Die Kenntnis der atomaren Struktur ist sehr wichtig um die physikahschen Eigen- 
schaften amorpher Materiahen zu verstehen. Die Strukturen ungeordneter Festkorper 
hangen normalerweise von ihren Praparationsmethoden ab. Daher is es wichtig, 
die Praparation zuerst zu modelheren. Dadurch konnen die resultierenden Struk- 
turen ermittelt werden. In Kapitel |^ wurden Niedrigenergie Molecular-Dynamic 
Simulationen von Atomstrahlwachstum von amorphem Kohlenstoff auf einer Dia- 
mant [111] Oberflache durchgefiihrt. Unser Ziel war es, Modelle von amorphem 
Kohlenstoff mit iiber 100 Atomen mit einer quantenmechanischen Behandlung des 
interatomaren Potentials zu erhalten. Als Wechselwirkung wurde ein wohlgetestetes 
und weitbenutztes tight-binding Potential verwendet. Die Strukturen wurden durch 
atomweise Anlagerung auf dem Substrat gewachsen. Der Wachstumsprozess wurde 
dem Experiment entsprechend ohne die Annahme einer kiinstlichen Energiedissipa- 
tion simuliert. Wir haben zwei unterschiedliche mittlere Atomstrahlenergien Efteam 
= 1 und 5 eV benutzt, und auch zwei veschiedene Substrattemperaturen Tsub = 
100 und 300 K fiir eine Simulationszeit von 30 ps. Die Simulationen wurden um 
15 ps bei Tsub = 100 K verlangert um grof^ere Strukturen herzustellen. Es wur- 
den sechs Netzwerke in zwei Dimensionen mit periodischen Randbedingungen durch 
diese Depositionstechnik erzeugt. AnschlieBend wurden die wichtigsten strukturellen 
Parameter, wie die Dichten, radialen Verteilungsfunktionen, Koordinationszahlen, 
die Verteilungsfunktionen der Bindungswinkel und Ringstatistiken analysiert. Die 
Depositionsgeschwindigkeit war groBer als in den Experimenten. Auf Grund un- 
serer Untersuchung zur zeitlichen Entwicklung sind wir sicher, dafi eine niedrigere 
Geschwindigkeit der Anlagerung keine wesentlichen Unterschiede in unseren Mod- 
ellen verursachen wiirde. Unsere niedrigeren Substrattemperaturen verglichen mit 
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den Experimenten stellen ein schnellere Energiedissipation sicher, die die hohe Depo- 
sitionsgeschwindingkeit einigermafien kompensiert. 

In Kapitel haben wir die atomistische Simulation des Chemical- Vapor-Deposi- 
tionsprozesses wahrend der Phase der bias-enhanced nucleation (BEN) von Dia- 
mant durchgefiihrt. Besonderes Augenmerk wurde auf die Untersuchung des Projek- 
tilbeschusses mit verschiedenen Hydrokarbonaten (CxHy) bei unterschiedlichen En- 
ergien gerichtet. Die Simulationen zeigten, daB groBere kinetische Energien tiefere 
Subplantationen zur Folge haben und daB es wahrscheinlicher wird, dafi diese Spezies 
dissoziert werden. Fiir Azetylen (C2H2) wird die Kohlenstoff-Kohlenstoff Bindung 
fast nie bei niedrigen Energien (20, 40 eV) gebrochen. Sie wird andererseits in 50 bis 
90% der Falle fiir groBere Energien (80 eV) gebrochen. Die Computer Experimente 
offenbarten, daB ungefahr 60% der Wasserstoff-Kohlenstoff Bindungen der Projektile 
wahrend der Subplantation gebrochen werden. Normaleweise wird bei Methyl (CH3) 
eine der Wasserstoff-Kohlenstoff Bindungen nicht gebrochen und die anderen beiden 
formen H2 Molekiile. Bei einer typischen bias-Spannung von 200 V, die optimal fiir 
die Formation von Diamantkeimen in BEN ist, haben die C2H2 Projektile 80 eV und 
die CH3 lonspezies 40 eV kinetische Energie. Es gibt bedeutensvoUe Unterschiede in 
den strukturellen Umordungen der Filme, die von den beiden Projektilen bei dieser 
Spannung verursacht wurden. Fiir die C2H2 Spezies ist die Durchdringung tiefer und 
diese Projektile sind verantwortlich fiir die Anderungen in tiefern Regionen als es der 
Fall mit den CH3 Spezies ist. Die CH3 Projektile tragen nicht zu Anderungen in tiefen 
Lagen der Filme bei und verursachen strukturelle Anderungen nur an der Oberflache. 
Die strukturellen Anderungen in den Filmen deuten darauf hin, daB Spezies mit hoher 
kinetischer Energie einen dickeren sp^-reichen Film unterhalb der mittleren Deposi- 
tionstiefe erzeugen und daB die Oberflache aus Atomen mit sp"^ Hybriden besteht. 
Die Grenze zwischen diesen beiden Regionen bewegt sich mit groBerer Depositionsen- 
ergie nach unten. Diese Ergebnisse unterstiitzen das Modell der Subplantation und 
deuten darauf hin, daB sp^ Hybride unter der Durchdringungschwelle geordnet werden 
konnen. Zukiinftig ware es eine sehr interessante Aufgabe, Simulationen zur Silizium- 
struktur zu unternehmen sowie die Auswirkungen dieser kristallinen Schicht auf die 
Formation von Diamantkeimen in den amorphen Kohlenstoffilmen zu untersuchen. 

Kapitel ^ handelt von der Simulation von amorphen Siliziumstrukturen. Genau 
das gleiche Programm fiir die Deposition, wie es in Kapitel Q fiir Kohlenstoff entwick- 
elt wurde, fand hierbei Anwendung. Wir haben einen gut getesteten tight-binding 
Hamiltonian benutzt um die interatomare Wechselwirkung zwischen Siliziumatomen 
zu beschreiben. Neun verschiedene amorphe Siliziumstrukturen wurden durch un- 
sere numerische Depositionstechnik generiert. Zusatzlich wurde ein zehntes Modell 
durch schnelle Abkiihlung prapariert, um den Unterschied zwischen atomarer Depo- 
sition und Abkiihlen aus der Schmelze zu untersuchen. Der wohl bedeutungsvoll- 
ste Unterschied zwischen diesen Modellen lag in den result ierenden Ringstatistiken, 
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beziehungsweise in der Medium- Range Order. Selbst Dreiecke sind in der atomaren 
Struktur sichtbar. Dies unterstiitzt drei andere unabhangige Hinweise auf die Exis- 
tenz von Dreiecken in amorphen Siliziumstrukturen. Dennoch wurden in der Literatur 
strukturelle Einheiten wie gleichseitige Dreiecke weder bei Rechnungen zur elektro- 
nischen Zustandsdichte noch bei der Theorie des Aufbrechens schwacher Bindungen 
beriicksichtigt. Letzterer ist der dominante Mechanismus bei der lichtinduzierten 
Defekterzeugung. Das tight-binding Modellpotential in unseren Computer Experi- 
menten iiberschatzt etwas die Bindungslangen, indem es Werte fiir Silizium in der 
fliissigen Phase hefert. Es erzeugt iiberkoordinierte Atome (fiinffach und sogar sechs- 
fach koordinierte Atome) wahrend der Simulation der amorphen Struktur. Die sechs- 
fach koordinierten Atome finden sich in der fliissigen Phase von Silizium, nicht jedoch 
in der amorphen Phase. Das wohlbekannte klassische, empirische Stillinger- Weber 
Potential fiir amorphes Silizium hatte das gleiche Problem und wurde kiirzlich mod- 
ifiziert. Es ware interessant, unser tight-binding Potential auf ahnliche Weise zu 
verbessern. 



In Kapitel ^ schlieBlich wurden Rechnungen fiir die Hopping-Driftbeweglichkeit 
von Ladungstragern in einer ein-dimensionalen Kette mit lokalisierten elektronischen 
Zustanden mit energetischer Unordnung ausgefiihrt. In den analytischen Rechnun- 
gen wurden exakte Ergebnisse wurden fiir Abhangigkeit der Beweglichkeit von der 
Temparatur und dem Feld fiir Systeme mit unkorrelierter Unordnung (Gaussches 
Unordnungsmodell) erzielt und soche fiir korrelierte Unordnug (korreliertes Unord- 
nungsmodell) fiir asymmetrische und symmetrische Ubergangswahrscheinlichkeiten. 
Es zeigte sich, daB diese Abhangigkeiten von der Lange des Systems abhangen. Daher 
konnten exakte Losungen, wie sie iiblicherweise fiir unendliche Systeme erzielt werden, 
fiir realistische diinne Systeme, wie sie in Anwendungen vorkommen, nicht relevant 
sein. In Bezug auf die Beschreibung der Transporteigenschaften von diskotischen 
fiiissigkristallinen Glasern sollte man daran erinnern, daB die echten Device-Proben 
normalerweise nur 50 bis 100 molekulare Sites in der Stapelrichtung aufweisen. Hi- 
erf iir sind die in Abschnitt ^]2| fiir endliche Systeme erzielten Ergebnisse relevant. In 
den analytischen Rechnungen haben wir nur Ubergange der Ladungstrager zu den 
Nachbarsites auf der Kette betrachtet. Die Resultate der Monte-Carlo Computersim- 
ulationen zeigen, daB dagegen fiir Systeme mit unkorrelierter raumlich-energetischer 
Unordnung Variable- Range- Hopping eine wichtige Rolle spielt. Langere Elektronen- 
spriinge konnen die Resultate fiir die Feld- und Temperaturabhangigkeiten der Be- 
weglichkeit der Ladungstrager andern. Fiir solche Systeme konnen die analytischen 
Ergebnisse fiir die Transportkoeffizienten (aus Teil |5]^) als eine Approximiation be- 
trachtet werden. Nimmt man zum Beispiel Spriinge iiber mehrere Nachbarn an, 
so bekommt man Werte fiir die Koeffizienten C, die die Temparaturabhangigkeit 
der Mobilitat beschreiben (|5.41|) , die ungefahr 10% von denjenigen abweichen, die 
man fiir ausschlieBliche Spriinge zu den nachsten Nachbarn erhalt. Mit einer solchen 
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Genauigkeit kann man sogar analytische Ergebnisse (aus Teil ^.2| ) fiir Systeme mit 
unkorrelierter Unordnung benutzen. Jedoch ist im Fall korrelierter raumlich-energe- 
tischer Unordnung dieser Effekt viel weniger ausgepragt. Dies impliziert, daB die 
exakten analytischen Ergebnisse fiir Hiipfen zum nachsten Nachbarn fiir solche Sys- 
teme ohne Einschrankung giiltig sind. 



Osszefoglalas 



Ez a dolgozat rendezetlen rendszerek kiilonbozo fizikai parametereinek szamitogepes 
modszerrekkel valo vizsgalataval foglalkozik. Az utobbi evtizedekben az informacios 
techno logiaban a hihetleniil inteziv es gyors fej lodes lehetove tette a fizikusok szamara 
a szamitogepek nagyon hatekony hasznalatat. A szeleskoru alkalmazasok kifejlesztet- 
tek egy lij es ma mar szeles korii reszet a fizikanak, amely valahol az elmeleti es 
ki'serleti fizika kozotti hatarmezsgyen talalhato. 

Az amorf anyagok fizikai tulajdonsagainak megertesehez nagyon fontos az atomi 
elrendezes ismerete. A rendezetlen rendszerek szerkezete altalaban fiigg az elkeszite- 
siik modszeretol. Ezert fontos, liogy ezt a fizikai folyamatot modelleziik eloszor, majd 
utana vizsgaljuk a kesz szerkezeteket. A |^. fejezetben, amorf szen [111] szubsztra- 
tum felett valo molekularis dinamikai novesztesenek kisenergias, atomi bombazasos 
modszer segitsegevel vizsgaltuk. Az volt a celunk, liogy tobb mint 100 reszecsket tar- 
talmazo, atomi skalan nagymeretii, amorf szen modelleket allitsunk elo a potencial 
kvantummechanikai leirasaval. A kolcsonhatast egy alaposan tesztelt es szeles korben 
alkalmazott szoros koteso (tight-binding) potenciallal modelleztiik. A szerkezetek az 
atomok a szubsztratumra valo egymas utani parologtatasos novesztesevel kesziiltek. 
A novekedes a kfserleteknek megfeleloen volt szimulalva. Ket kiilonbozo atlagos bom- 
bazasi energiat, Ef,eam = 1 es 5 eV, es ket kiilonbozo szubsztrat homersekletet Tsub = 
100 es 300 K alkalmaztunk a 30 ps hosszii szimulacios ido alatt. A szimulaciok 15 ps- 
mal meg voltak hosszabitva Tsub = 100 K eseten, hogy nagyobb meretii szerkezeteket 
kapjunk. Hat darab modell kesziilt el ezzel a parologtatasos modszerrel, ket— dimen- 
zios periodikus hatarfeltetelt alkalmazva. A legfontosabb szerkezeti parametereket, 
mint a siiriiseget, radialis eloszlas fiiggvenyt, koordinacios szamot, szogeloszlast, es 
gyiiriistatisztikat elemeztiik. A parologtatasi rata sokkal nagyobb volt mint a kiser- 
letekben. A rendszer idofejlodeset koveto vizsgalatunk alapjan azonban azt varjuk, 
hogy sokkal kisebb parologtatasi rata eseten sem kapnank jelentos valtozasokat a mo- 
delljeinkben. Az altalunk hasznalt, kiserletekhez kepest sokkal alacsonyabb szubszt- 
ratum homerseklet gyorsabb energia disszipaciot eredmenyez, es kisse kompenzalja a 
nagy parologtatasi ratat. 

A|^. fejezetben, gyemant, kemiai gozfazis levalasztasii (Chemical Vapor Depositi- 
on - CVD), elektromos elofesziteses nukleacios fazisanak (Bias Enhanced Nucleation - 
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BEN) atomi szintu modellezese van osszefoglalva. Kiilonos figyelmet forditottunk kii- 
lonbozo szenhidrat (C^H^) gyokok, mas-mas kinetikus energiaval valo bombazasara. 
Az eredmenyek azt mutatjak, hogy nagyobb kinetikus energia melyebb szubplantaci- 
6t okoz es a gyokok disszociacioja is gyakoribb. Az acetilen (C2H2) eseten a szen-szen 
kotes nagyon ritkan bomlik fel alacsony energas (20, 40 eV) bombazas eseten, viszont 
az esetek 50— 90%-ban felbomlik nagyobb energiak eseten (80 eV). A szami'togepes 
ki'serletek kimutattak, hogy koriilbeliil az esetek 60%-ban a szen-hidrogen kotes fel- 
bomlik a szubplantacio soran. A metil (CH3) eseten altalaban az egyik szen-hidrogen 
kotes megmarad, es a masik ket hidrogen H2 molekulat alkot. A tipikus 200 V-os elo- 
feszites eseten, ami optimalis a gyemant nukleusz kepzodesehez a BEN soran, a C2H2 
gyokok 80 eV-os, mig a CH3 ion gyokok 40 eV-os kinetikus energiaval rendelkeznek. 
Ennel az elofeszitesnel, jelentos kiilonbsegek figyelhetoek meg a ket gyok altal oko- 
zott amorf szen vekonyretegbeli szerkezeti atrendezodeseben. A C2H2 gyok eseten, 
a behatolas melyebb es a gyok melyebb tartomanyokban okoz valtozasokat mint a 
CH3 gyok. A CH3 lovedek nem jarul hozza mely valtozasokhoz az amorf retegben es 
csak a feliilet kozeleben okoz struktiiralis atalakulasokat. A vekonyretegben megfi- 
gyelt szerkezeti valtozasok kimutattak, hogy nagyobb kinetikus energiajii gyokok az 
atlagos behatolasi melyseg alatt egy vastag, sp^ hibridekben gazdag reteget hoznak 
letre es a feliilet pedig sp^ hibridekkel rendelkezo atomokbol all. A hatar ekozott a 
ket reteg kozott lejjebb tolodik nagyobb bombazasi energia eseten. Az eredmenye- 
ink megerosftik a szubplantacios modell alkalmazasat. Egy igen erdekes feladat lehet 
a jovoben megismetelni a szimulacios ki'serleteket egy szilfcium szubsztratum felett, 
es vizsgalni ennek a kristalyos retegnek a hatasat az amorf szen vekonyretegben a 
gyemant nukleuszok kialakulasahoz vezeto folyamatat. 

A ^ fejezet az amorf szilfcium szerkezetenek szimulaciojaval foglalkozik. Ugya- 
nazt az atomi parologtatasi programcsomagot hasznaltuk, mint amelyet kifejlesztet- 
tiink az amorf szen novesztesenek es a 0. fejezetben mutattunk be. 

Egy szeles korben tesztelt, szoros kotesii (tight-binding) Hamiltont formalizmust al- 
kalmaztunk a szilfcium atomok kozotti kolcsonhatas lefrasara. Egy hasonlo szoros 
koteso potencialt fejlesztett ki ugyanez a kutatocsoport a szen rendszerekre es mi ezt 
alkalmaztuk eredmenyesen az amorf szen vekonyretegek novesztesenek lefrasara (Id. 
^ fejezet). Kilenc kiilonbozo amorf szilfcium szerkezetet keszftettiink el a szamfto- 
gepes parologtatasi modszeriinkkel. Egy tovabbi tizedik modell gyorshiiteses techni- 
kaval kesziilt el, hogy osszehasonlftsuk az atomi parologtatasi es az olvadek hfiteses 
modszerek kozotti kiilonbsegeket. Az altalunk talalt legfontosabb kiilonbseg a gyfi- 
riistatisztikaban figyelheto meg, azaz a kozeptavu rendben. Haromszogek ugyancsak 
megtalalhatoak az atomi elrendezesben, megerosftve ezzel harom korabbi fiiggetlen 
bizonyossagat a haromszogek elofordulasanak az amorf szilfciumban. Mindamellett, 
olyan fragmensek, mint peldaul egyenlo oldalu haromszogek sohasem voltak figyelem- 
be veve az elektron allapotsfirfiseg szamftasanal, illetve olyan megoldatlan problemak 
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eseten, mint a gyenge kotesek felbomlasa hosszii megvilagitas eseten, ami a fo me- 
chanizmusa a feny altal indukalt hiba keltesenek. A szoros kotesii modell potencial 
egy picit tiilbecsiili a kotestavolsagot a szami'togepes ki'serleteinkben, olyan ertekeket 
adva, amelyeket korabban folyadek sziliciumban figyeltek meg. Az amorf szilfcium 
szerkezetenek szimulacioja soran a potencial tiilkoordinalt atomokat is ad (otos, sot 
hatos koordinaltsagu atomokat is). Hatos koordinaltsagu atomok, inkabb a szilfci- 
um folyadek fazisaban talalhatoak, mint az amorf sziliciumban. Az amorf szilfcium 
eseteben a jol ismert Stillinger- Weber klasszikus empirikus potencial is liasonlo ne- 
hezsegekkel kiizdott, es mostanaban eredmenyesen lett modosftva. Egy jovobeli eel 
lehet a szoros kotesfi (tight-binding) potencial hasonlo modon valo javftasa. 

Vegezetiil, az ||. fejezetben, egy egy-dimenzios rendezetlen lokalizacios energia 
allapotokon keresztiil lialado tolteshordozok mozgekonysaganak szamftasa talalhato 
meg. Az analitikus szamftasokban, egzakt eredmenyeket kaptunk a mozgekonysag 
elektromos tererosseg es homerseklet fiiggesere mind nemkorrelalt rendezetlen rend- 
szer eseten (Gaussian disordered model - GDM) es mind korrelacio eseten is (corre- 
lated disorder model - GDM), asszimetrikus es szimmetrikus atmeneti ratak eseten 
egyarant. Ezek az eredmenyek a rendszer meretetol fiiggenek. Ezert meg az eg- 
zakt megoldasok vegtelen hosszii rendszerek eseten is ervenyiiket veszthetik olyan 
vekony mintak eseten, amelyeket a valodi alkalmazasokban hasznalnak. Kozeppont- 
ba helyezve az ugy nevezett oszloposan korongszerfi folyadekkristalyokat, figyelembe 
kell venni, hogy az alkalmazott eszkozok altalaban csak 50-100 "epfto elembol" all- 
nak. Ilyen rendszerek eseten az f>.2\ alfejezetben kapott elmeleti eredmenyek koziil 
a veges meretfi rendszerekre kapott eredmenyeket kell alkalmazni, a vegtelen hosszu 
rendszerekre levezetett eredmenyek helyett. Az analitikus szamftasok soran a toltes- 
hordozok hopping] anal csak elsoszomszed atmeneteket vettiink figyelembe. A Monte 
Garlo szimulaciok eredmenyei kimutattak, hogy azokban a rendszerekben, ahol nines 
korrelacio az energia es lokalizacios allapotok terbeli elhelyezkedese kozott, a valtozo 
tavolsagii hopping fontos szerepet jatszik es hatassal van a tolteshordozo mozgekony- 
saganak tererosseg es homerseklet fiiggesre (Id. |5.3| . alfejezet). Ezekre a rendszerekre 
a transzport egyiitthatokra kapott analitikus eredmenyek kozelfteskent alkalmazha- 
toak. Peldaul, ha tavolabbi szomszed hoppingot is figyelembe vesziik, akkor a C 
egyiitthatora, ami a mozgekonysag homerseklet fiiggeset frja le az ( p.41| ) egyenletben, 
10%-os kiilonbseg talalhato a csak elsoszomszed hoppingra levezetett eredmenyekhez 
kepest. Ilyen hiba mellett, akar az analitikus eredmenyek (Id. |5.2| . alfejezet) is merva- 
doak lehetnek a korrelalatlan rendszerekben. Ugyanakkor, a korrelalt rendezetlenseg 
eseten, ez az effektus sokkal kisebb, utalva arra, hogy az elsoszomszed hoppingra ka- 
pott egzakt analitikus eredmenyek kiilonosebb korlatozasok nelkiil ervenyesek ezekben 
a rendszerekben. 
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